IPPP/22/71, DCPT/22/142

# High Energy Resummed Predictions for the Production of a Higgs Boson with at
least One Jet

Jeppe R. Andersena, Hitham Hassana, Andreas Maierb, Jérémy Paltrinieric,
Andreas Papaefstathioud and Jennifer M. Smilliec

a Institute for Particle Physics Phenomenology,
University of Durham, South Road, Durham DH1 3LE, UK
b Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany
c Higgs Centre for Theoretical Physics, University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh EH9 3FD, UK,
dDepartment of Physics, Kennesaw State University, Kennesaw, GA 30144, USA

###### Abstract

We present all-order predictions for Higgs boson production plus at least one
jet which are accurate to leading logarithm in $\hat{s}/|p_{\perp}|^{2}$. Our
calculation includes full top and bottom quark mass dependence at all orders
in the logarithmic part, and to highest available order in the tree-level
matching. The calculation is implemented in the framework of High Energy Jets
(HEJ). This is the first cross section calculated with $\log(\hat{s})$
resummation and matched to fixed order for a process requiring just one jet,
and our results also extend the region of resummation for processes with two
jets or more. This is possible because the resummation is performed explicitly
in phase space. We compare the results of our new calculation to LHC data and
to next-to-leading order predictions and find a numerically significant impact
of the logarithmic corrections in the shape of key distributions, which
remains after normalisation of the cross section.

###### Contents

  1. 1 Introduction
  2. 2 Higgs Boson plus Jets Production in the High-Energy Limit
    1. 2.1 Scaling of Amplitudes at High Energies
    2. 2.2 $H+\geq 2j$ Processes within HEJ
    3. 2.3 Scaling of $H+\geq 1j$ Amplitudes
    4. 2.4 New Components for $H+\geq 1j$ and an Outer Higgs Boson
      1. 2.4.1 Higgs Current
      2. 2.4.2 Lipatov Vertex for Additional Gluons
      3. 2.4.3 Matrix element including additional gluons
  3. 3 Predictions and Comparison to Data
    1. 3.1 Predictions
    2. 3.2 Predictions for 13 TeV and Comparison to Data
    3. 3.3 Predictions for 8 TeV and Comparison to Data
      1. 3.3.1 $H+\geq 1j$
      2. 3.3.2 $H+\geq 2j$
  4. 4 Conclusions
  5. 5 Appendix - NLO reweighting factors

##  1 Introduction

Analysing the Higgs sector is among the foremost objectives of the LHC. To
this end, experiments aim for accurate measurements of processes where Higgs
bosons are produced, either inclusively or in association with other
identified particles. Given the phenomenological importance of processes
involving Higgs boson production, there are considerable efforts to provide
high-precision theory predictions. Perturbative corrections are typically
found to be sizable, necessitating the inclusion of effects at high orders.
This endeavour faces a major challenge: in large regions of phase space Higgs
boson production is predominantly loop-induced, namely through gluon fusion
via a virtual top-quark loop. Inclusive gluon-fusion Higgs boson production
with full finite top-mass contributions is currently known at next-to-next-to-
leading order (NNLO) [1], exclusive Higgs boson plus jet production at next-
to-leading order (NLO) [2, 3], and Higgs boson plus dijet production only at
leading order (LO) [4, 5].

To facilitate calculations, the top-quark mass is often assumed to be much
larger than all other scales. Based on this approximation, one more order has
been computed in the perturbative expansion for the aforementioned processes
[6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. However, one is often interested in
observables where the assumption of a comparatively large top-quark mass is
invalid and the full mass dependence has to be accounted for. One example is
the study of the high-energy tail in the Higgs boson transverse momentum
distribution.

Another avenue towards better theory predictions consists of the all-order
resummation of contributions that are enhanced in kinematic regions of
interest. For Higgs boson production together with at least one jet, one finds
logarithms in $\hat{s}/|p_{\perp}|^{2}$, where $\hat{s}$ is the square of the
partonic centre-of-mass energy and $p_{\perp}$ a characteristic transverse
momentum scale [16]. For the case of two or more jets, the resummation of
these high-energy logarithms has been shown to lead to significant
corrections, especially after weak-boson fusion cuts are applied [17]. This
provides a strong motivation to extend the study of logarithmic enhancement to
the production of a Higgs boson with a single jet.

Inclusive calculations of BFKL resummation for Higgs boson plus jet production
have been performed [18, 19]. In contrast, our resummation of high-energy
logarithms is based on the _High Energy Jets_ (HEJ) framework [20, 21, 22,
23]. HEJ provides realistic predictions through a fully flexible Monte Carlo
implementation, supplementing leading-order perturbation theory with high-
energy resummation retaining exact gauge invariance and momentum conservation.
The calculation presented here is the first time this approach has been used
for an inclusive 1-jet process. As is necessary in the high-energy region, the
all-order resummation includes the full effects of finite quark masses. We
first review the formalism and derive the new building blocks required for
leading-logarithmic (LL) resummation for Higgs boson plus jet production in
section 2. In section 3, we compare our predictions to experimental
measurements and propose observables tailored to the systematic analysis of
high-energy corrections. We conclude in section 4.

##  2 Higgs Boson plus Jets Production in the High-Energy Limit

In the following, we discuss the general properties and structure of
amplitudes in the high-energy limit. We briefly summarise LL resummation in
the _High Energy Jets_ formalism and derive the new ingredients for the
production of a Higgs boson together with a single jet, and for processes with
two or more jets where the Higgs boson is outside of the jets.

###  2.1 Scaling of Amplitudes at High Energies

Generally, we are interested in the behaviour of amplitudes in the region of
_Multi-Regge Kinematics_ (MRK). This region is defined by a large centre-of-
mass energy with large invariant masses between all pairs of outgoing
particles with finite transverse momenta. This is equivalent to a strong
ordering in rapidities. Specifically, for a $2\rightarrow n$ process, we
require

$\displaystyle y(p_{n})\gg...\gg
y(p_{1})\qquad|p_{i\perp}|\sim\text{finite}\;\forall i\in\\{1..n\\},$ (1)

where the outgoing particle $i$ has momentum $p_{i}$, rapidity $y_{i}\equiv
y(p_{i})$ and transverse momentum $|p_{i\perp}|$.

In this region, Regge theory [24] states that the amplitude should scale as

$\mathcal{M}\sim s_{12}^{\alpha_{1}(t_{1})}...s_{n\,n+1}^{\alpha_{n}(t_{n})},$
(2)

where the $s_{i\,i+1}$ refers to the invariant mass between particle $i$ and
$i+1$, and $\alpha_{i}(t_{i})$ is the maximum spin of any particle that can be
exchanged in the $t$-channel between particle $i$ and $i+1$. From this
formula, it follows that the leading contribution to a QCD amplitude is given
by the configurations which maximise the number of gluons exchanged in the
$t$-channel. These configurations characterise the regions of phase space in
which the leading high-energy logarithms arise. We therefore refer to them as
leading-logarithmic (LL) or Fadin-Kuarev-Lipatov (FKL) configurations.

As a simple example, let us consider the amplitude for elastic scattering of a
quark or antiquark ($q$) and a gluon ($g$), with the incoming quark in the
backward direction [21]. Ordering the outgoing particles by ascending
rapidity, the two possible configurations are $qg\to qg$ and $qg\to gq$. For
the rapidity ordering $qg\to qg$ it is possible to exchange a $t$-channel
gluon and we therefore expect the amplitude to scale as $\mathcal{M}_{qg\to
qg}\sim\hat{s}^{1}$ for $y(p_{2})\gg y(p_{1})$. Conversely, the flipped
ordering $qg\to gq$ only allows a $t$-channel (anti-)quark exchange, implying
$\mathcal{M}_{qg\to gq}\sim\hat{s}^{\frac{1}{2}}$. This scaling behaviour is
indeed confirmed by an explicit calculation and illustrated in figure 1, where
increasing $\Delta y$ represents approaching the MRK limit.

Figure 1: Regge scaling of quark-gluon scattering amplitudes in the MRK limit,
approached by increasing $\Delta y$. Left: Rapidity ordering $qg\to qg$.
Right: Rapidity ordering $qg\to gq$. The kinematics are fixed by the azimuthal
angle $\phi_{1}=\frac{\pi}{7}$ and transverse momentum $p_{1\perp}=40\,$GeV of
the outgoing particle in the backward direction.

###  2.2 $H+\geq 2j$ Processes within HEJ

The construction of the leading-logarithmic calculation of $pp\to H+\geq 2j$
in the HEJ framework was described in detail in [25, 17]. Here we summarise
the main points in order to frame the discussion of the new components
calculated in this paper.

Following the arguments in section 2.1, the LL configurations in pure QCD have
the form $f_{a}f_{b}\to f_{a}\cdots f_{b}$, where $f_{a},f_{b}$ indicate the
incoming parton flavours and the ellipsis denotes an arbitrary number of
gluons. As before, the particles are written in order of increasing rapidity.
The production of an additional Higgs boson proceeds via an effective coupling
to two or more gluons. Since invariant masses are large in the high-energy
region, it is crucial that the exact dependence on the top-quark mass is
included in this effective coupling. A final-state Higgs boson with momentum
$p_{H}\equiv p_{j}$ at an intermediate rapidity $y_{j}$ such that $y_{j-1}\ll
y_{j}\ll y_{j+1}$ can then exchange $t$-channel gluons with the outgoing
partons $j-1,j+1$. It was shown in [25] that the scaling behaviour in equation
(2) directly generalises when a Higgs boson is emitted in the middle of the
quarks and gluons. Therefore, all configurations $f_{a}f_{b}\to f_{a}\cdots
H\cdots f_{b}$ contribute at LL accuracy.

In the MRK limit the amplitudes are found to factorise into a neat product of
simple functions. In the _High Energy Jets_ formalism we obtain the form

$\begin{split}\overline{\left|\mathcal{M}_{{\tt HEJ}}^{f_{a}f_{b}\to
f_{a}\cdots H\cdots
f_{b}}\right|}^{2}={}&\mathcal{B}_{f_{a},H,f_{b}}(p_{a},p_{b},p_{1},p_{n},q_{j},q_{j+1})\\\
&\cdot\prod_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{n-2}\mathcal{V}(p_{a},p_{b},p_{1},p_{n},q_{i},q_{i+1})\\\
&\cdot\prod_{i=1}^{n-1}\mathcal{W}(q_{i},y_{i},y_{i+1}),\end{split}$ (3)

for the modulus square of the matrix element, summed and averaged over
helicities and colours. In this expression, $p_{a}$ $(p_{b})$ is the incoming
momentum in the backward (forward) direction and $p_{1},\dots,p_{n}$ are the
outgoing momenta ordered in increasing rapidity. The $t$-channel momenta are
given by

$q_{1}=p_{1}-p_{a},\qquad q_{i}=q_{i-1}-p_{i}\;\text{ for }i>1.$ (4)

The structure is illustrated in figure 2.

Figure 2: Structure of the matrix element for the process $f_{a}f_{b}\to
f_{a}\cdots H\cdots f_{b}$.

At Born level, the right-hand side of equation (3) reduces to the function
$\mathcal{B}_{f_{a},H,f_{b}}$, described below. $\mathcal{V}$ comprises the
real corrections due to the production of $n-3$ gluons in addition to
$f_{a},f_{b}$, and the Higgs boson. It is given by the contraction of two
Lipatov vertices [25]:

$\displaystyle\mathcal{V}(p_{a},p_{b},p_{1},p_{n},q_{i},q_{i+1})={}$
$\displaystyle-\frac{C_{A}}{t_{i}t_{i+1}}V_{\mu}(p_{a},p_{b},p_{1},p_{n},q_{i},q_{i+1})V^{\mu}(p_{a},p_{b},p_{1},p_{n},q_{i},q_{i+1}),$
(5) $\displaystyle V^{\mu}(p_{a},p_{b},p_{1},p_{n},q_{i},q_{i+1})={}$
$\displaystyle-(q_{i}+q_{i+1})^{\mu}$
$\displaystyle+\frac{p_{a}^{\mu}}{2}\left(\frac{q_{i}^{2}}{p_{i+1}\cdot
p_{a}}+\frac{p_{i+1}\cdot p_{b}}{p_{a}\cdot p_{b}}+\frac{p_{i+1}\cdot
p_{n}}{p_{a}\cdot p_{n}}\right)+p_{a}\leftrightarrow p_{1}$
$\displaystyle-\frac{p_{b}^{\mu}}{2}\left(\frac{q_{i+1}^{2}}{p_{i+1}\cdot
p_{b}}+\frac{p_{i+1}\cdot p_{a}}{p_{b}\cdot p_{a}}+\frac{p_{i+1}\cdot
p_{1}}{p_{b}\cdot p_{1}}\right)-p_{b}\leftrightarrow p_{n},$ (6)

where $t_{i}=q_{i}^{2}$ are the squares of the $t$-channel momenta.
$\mathcal{W}$ accounts for the all-order finite contribution coming from the
sum of the virtual corrections and unresolved real corrections. It is process-
independent and described in detail in [25].

The process-dependent Born-level factor is given by

$\mathcal{B}_{f_{a},H,f_{b}}=\frac{(4\pi\alpha_{s})^{n-1}}{4(N_{c}^{2}-1)}\frac{K_{f_{a}}(p_{1}^{-},p_{a}^{-})}{t_{1}}\
\frac{K_{f_{b}}(p_{n}^{+},p_{b}^{+})}{t_{n-1}}\frac{\left\|S_{f_{a}f_{b}\to
f_{a}Hf_{b}}\right\|^{2}}{t_{j}t_{j+1}}.$ (7)

Here, $\alpha_{s}$ is the strong coupling and $N_{c}=3$ is the number of
colours. The difference between incoming gluons and (anti-)quarks is
completely absorbed into the _colour acceleration multipliers_ $K_{f}$ with

$\displaystyle K_{g}(x,y)={}$
$\displaystyle\frac{1}{2}\left(\frac{x}{y}+\frac{y}{x}\right)\left(C_{A}-\frac{1}{C_{A}}\right)+\frac{1}{C_{A}}$
for gluons, (8) $\displaystyle K_{q}(x,y)={}$ $\displaystyle C_{F}$ for quarks
and antiquarks. (9)

$C_{F}=\frac{N_{C}^{2}-1}{2N_{C}}$ and $C_{A}=N_{C}$ are the usual Casimir
invariants. $S_{f_{a}f_{b}\to f_{a}Hf_{b}}$ is a contraction of currents with
the Higgs boson production vertex. The double vertical bars indicate the sum
over helicities of the corresponding amplitudes:

$\left\|S_{f_{a}f_{b}\to
f_{a}Hf_{b}}\right\|^{2}=\sum_{\begin{subarray}{c}\lambda_{a}=+,-\\\
\lambda_{b}=+,-\end{subarray}}\left|j^{\lambda_{a}}_{\mu}(p_{1},p_{a})V_{H}^{\mu\nu}(q_{j},q_{j+1})j_{\nu}^{\lambda_{b}}(p_{n},p_{b})\right|^{2}.$
(10)

$V_{H}$ is the well-known one-loop effective coupling between the Higgs boson
and two gluons in the normalisation of [17], including the full quark-mass
dependence. The inclusion of this piece in equation (3) then gives the correct
finite quark-mass contributions at LL for _any_ number of final state
partons/jets. Finally, the current $j$ is given by

$j^{\lambda}_{\mu}(p,q)=\bar{u}^{\lambda}(p)\gamma_{\mu}u^{\lambda}(q).$ (11)

In addition to the LL resummation discussed so far, gauge-invariant subsets of
next-to-leading logarithmic (NLL) corrections originating from non-FKL
configurations have also been included in HEJ. One source of NLL corrections
are the configurations $qf_{b}\to Hq\cdots f_{b}$ and $f_{a}q\to f_{a}\cdots
qH$, which only permit $n-2$ $t$-channel gluon exchanges instead of the $n-1$
exchanges found in LL configurations. In these cases, we adapt the matrix
element formula for the corresponding LL configurations to a flipped rapidity
order of outgoing (anti-)quark and Higgs boson. If the Higgs boson is emitted
first in rapidity order, we use equation (3) with $p_{2}\equiv p_{H}$ and
exclude the virtual correction factor $\mathcal{W}$ for $i=1$. In the other
case of the Higgs boson being emitted last, we set $p_{n-1}=p_{H}$ and skip
$\mathcal{W}$ for $i=n-1$.

A second class of non-FKL configurations arises for three or more produced
jets, when the most backward or forward outgoing particle is a gluon, but the
corresponding incoming parton is a quark or antiquark. These “unordered gluon”
configurations, $qf_{b}\to gq\cdots H\cdots f_{b}$ and $f_{a}q\to f_{a}\cdots
H\cdots qg$, allow one $t$-channel gluon exchange less than the corresponding
FKL configurations in which the unordered gluon is swapped with the
neighbouring (anti-)quark. Hence, they contribute at NLL accuracy. Without
loss of generality, we consider the case where the unordered gluon is the most
backward emitted particle. We denote its momentum by $p_{g}$ and the following
momenta by $p_{1},\dots,p_{n}$. The modulus square of the matrix element then
has the same structure as in equation (3). In fact, the only changes are that
the first $t$-channel momentum is now $q_{1}=p_{a}-p_{1}-p_{g}$ and that a
different Born-level function $\mathcal{B}_{gq,H,f_{b}}$ depending also on
$p_{g}$ appears. For a derivation and explicit expressions, see [25].

###  2.3 Scaling of $H+\geq 1j$ Amplitudes

To extend the formalism to the production of a Higgs boson with a single jet
we first need to identify the LL configurations, following the discussion in
section 2.1, and then derive the corresponding matrix elements.

So far, we have only considered LL configurations in which both the most
backward and the most forward outgoing particle is a parton. However, in the
process $gq\rightarrow Hq$, the amplitude should scale as $\mathcal{M}\sim
s_{Hq}$, as there is a gluon exchange (thus a spin-1 particle) in the
$t$-channel. Similarly, the process $gg\rightarrow Hg$ corresponds to
$\mathcal{M}\sim s_{Hg}$. If we look at Higgs boson plus dijet production, the
same argument allows us to establish that $gq\rightarrow Hgq$ scales as
$\mathcal{M}\sim s_{Hg}s_{gq}$. All these configurations therefore contribute
at LL accuracy. This is no longer the case if, for example, outgoing parton
flavours are rearranged: $gq\rightarrow Hqg$ scales as $\mathcal{M}\sim
s_{Hq}s_{qg}^{1/2}$.

Note that these scalings are valid whether we consider the full LO amplitude
(with Higgs to gluons couplings _via_ quark loops) or the Higgs Effective
Field Theory (HEFT) one with an infinite top mass $m_{t}$, as shown in figure
3. To produce these plots, the amplitude is extracted from Madgraph5_aMC@NLO
[26] and is calculated in a one-dimensional phase-space as a function of the
rapidity separation between all pairs of particles. It was checked for
consistency that setting the top mass to infinity in the LO amplitude yields
the HEFT result. We compare to the LO truncation of the all-order HEJ
amplitudes, anticipating their derivation from the high-energy limit in
section 2.4,

(a) $|\mathcal{M}|^{2}/\hat{s}^{2}$ for the process $gu\rightarrow Hu$

(b) $|\mathcal{M}|^{2}/\hat{s}^{2}$ for the process $gg\rightarrow Hg$

(c) $|\mathcal{M}|^{2}/(s_{Hg}^{2}\times s_{gu}^{2}))$ for the process
$gu\rightarrow Hgu$

(d) $|\mathcal{M}|^{2}/(s_{Hu}^{2}\times s_{ug}))$ for the process
$gu\rightarrow Hug$

Figure 3: Verifying Regge scaling of the squared matrix elements (equation
(2)) for 4 different processes. Increasing values of $\Delta y$ on the
$x$-axis approach the MRK limit (equation (1)).

The momentum configurations chosen are summarised in table 1. We stress though
that the behaviour shown is not dependent on specific values of azimuthal
angle or transverse momentum, but only on the rapidity assignment of the
particles.

Process | Momenta configuration
---|---
$gq\rightarrow Hq$ | $\begin{cases}y_{q}=-\Delta\text{ and }y_{H}=\Delta&\\\ \phi_{q}=\frac{\pi}{4}&\\\ p_{q\perp}=40\,\text{GeV}&\end{cases}$
$gg\rightarrow Hg$ | $\begin{cases}y_{g}=-\Delta\text{ and }y_{H}=\Delta&\\\ \phi_{g}=\frac{\pi}{4}&\\\ p_{g\perp}=40\,\text{GeV}&\end{cases}$
$gq\rightarrow Hgq$ | $\begin{cases}y_{q}=-\Delta,y_{g}=0\text{ and }y_{H}=\Delta&\\\ \phi_{g}=\frac{\pi}{2}\text{ and }\phi_{q}=-\frac{\pi}{3}&\\\ p_{g\perp}=k_{q\perp}=40\,\text{GeV}&\end{cases}$
$gq\rightarrow Hqg$ | $\begin{cases}y_{g}=-\Delta,y_{q}=0\text{ and }y_{H}=\Delta&\\\ \phi_{g}=-\frac{\pi}{2}\text{ and }\phi_{q}=\frac{\pi}{3}&\\\ p_{g\perp}=k_{q\perp}=40\,\text{GeV}&\end{cases}$
Table 1: The momentum configurations used in figure 3.

###  2.4 New Components for $H+\geq 1j$ and an Outer Higgs Boson

In section 2.1 we discussed the factorisation of LL amplitudes for
$f_{a}f_{b}\to f_{a}\cdots H\cdots f_{b}$ into a Born-level function
$\mathcal{B}$, a product over real-emission vertices $\mathcal{V}$, and a
product of virtual corrections $\mathcal{W}$. The same type of factorisation
holds for LL configurations with the Higgs boson as the most forward or
backward outgoing particle. In fact, the virtual corrections are the same as
in equation (3). To derive the remaining factors, we first analyse the Born-
level process $gf_{b}\to Hf_{b}$ and then consider real corrections.

####  2.4.1 Higgs Current

The Born-level function $\mathcal{B}_{H,f_{b}}$ for the process $gf_{b}\to
H\cdots f_{b}$ is obtained by deriving a $t$-channel factorised form analogous
to equation (7) from the modulus square of the Born-level amplitude in the MRK
limit. For $gq\to Hq$, the tree-level amplitude is determined by a single
diagram, depicted in figure 4.

Figure 4: The tree-level diagram for the process $gq\to Hq$. The straight
solid line denotes an arbitrary light quark or antiquark.

Without requiring any approximations we obtain the factorised expression

$\displaystyle\mathcal{B}_{H,f_{b}}={}$
$\displaystyle\frac{(4\pi\alpha_{s})^{n-1}}{4(N_{c}^{2}-1)}\frac{1}{t_{1}}\frac{K_{f_{b}}(p_{n}^{+},p_{b}^{+})}{t_{n-1}}\left\|S_{gf_{b}\to
Hf_{b}}\right\|^{2},$ (12) $\displaystyle\left\|S_{gf_{b}\to
Hf_{b}}\right\|^{2}={}$
$\displaystyle\sum_{\begin{subarray}{c}\lambda_{a}=+,-\\\
\lambda_{b}=+,-\end{subarray}}\left|\epsilon_{\mu}^{\lambda_{a}}(p_{a})\
V_{H}^{\mu\nu}(p_{a},p_{a}-p_{1})\
j_{\nu}^{\lambda_{b}}(p_{n},p_{b})\right|^{2},$ (13)

for $f_{b}=q$, where $\epsilon^{\lambda_{a}}(p_{a})$ is the polarisation
vector of the incoming gluon. This is plotted along with the exact LO results
from Madgraph5_aMC@NLO [26] in figure 3(a), showing exact agreement for both
finite top quark mass and in the infinite $m_{t}$ limit. In the MRK limit,
this formula also holds for $f_{b}=g$, which is shown in figure 3(b). In this
case there is some approximation away from the limit, but very quickly the LO
and HEJ lines converge as $\Delta y$ increases.

####  2.4.2 Lipatov Vertex for Additional Gluons

In section 2.2, we described the simple factorised structure of amplitudes
within (N)MRK limits. Not only are the different components independent of
momenta in different parts of the chain, they are independent of the particle
content of the rest of the chain. This should mean that the Lipatov vertex
derived in pure QCD processes for additional gluons still applies. However,
the Lorentz and colour structure of the “Higgs current”
$j_{H}^{\nu}=\epsilon_{\mu}V_{H}^{\mu\nu}$ differ compared to pure QCD
processes so it is important to check that this is indeed the case.

|  |  |
---|---|---|---
a) | b) | c) | d)

|  |  |
---|---|---|---
e) | f) | g) | h)

Figure 5: The eight LO diagrams which contribute to the process in equation
(14).

We will consider the process

$\displaystyle g(p_{a})q(p_{b})\rightarrow H(p_{1})g(p_{2})q(p_{3}),$ (14)

in the MRK limit $y_{1}\ll y_{2}\ll y_{3}$. There are eight LO diagrams, as
shown in fig. 5. Compact expressions for tree-level Higgs-plus-4 parton
colour-ordered amplitudes appear in [27, 28]. Setting $q_{1}=p_{a}-p_{1}$ and
$q_{2}=p_{3}-p_{b}$, the HEJ amplitude is given by

$\displaystyle\begin{split}\mathcal{M}_{{\tt
HEJ}}=ig_{s}^{2}f^{2ea}T^{e}_{3b}\
\frac{\bar{u}(p_{3})\gamma^{\nu}u(p_{b})}{q_{1}^{2}q_{2}^{2}}\
V^{\alpha}(p_{a},p_{b},p_{a},p_{3},q_{1},q_{2})V_{H}^{\mu\nu}(p_{a},q_{1})\
\varepsilon_{\mu}(p_{a})\varepsilon_{\alpha}(p_{2})^{*}.\end{split}$ (15)

As the outer particle is no longer colour-charged, the third argument of the
Lipatov vertex defined in equation (2.2) is now $p_{a}$ instead of $p_{1}$.
The colour factor of the HEJ amplitude may be rewritten

$\displaystyle if^{2ea}T^{e}_{3b}=(T^{a}T^{2})_{3b}-(T^{2}T^{a})_{3b}.$ (16)

We can then directly compare equation (15) with the MRK limit of eqs. (26) and
(27) in Ref. [28], and we find agreement at LL up to an unphysical phase
arising from our spinor conventions. Specifically, the LL term in the MRK and
infinite top-quark mass limit of equation (15) is given by

$\displaystyle g_{s}^{2}\left(\frac{\alpha_{s}}{3\pi
v}\right)if^{2ea}T^{e}_{3b}\frac{\langle
3a\rangle[ab]}{|p_{2\perp}||p_{3\perp}|}\rightarrow
g_{s}^{2}\left(\frac{\alpha_{s}}{3\pi
v}\right)if^{2ea}T^{e}_{3b}\frac{\hat{s}}{|p_{2\perp}||p_{3\perp}|},$ (17)

where the angle and square brackets are Lorentz-invariant kinematic factors
defined by $\langle ij\rangle=\bar{u}(p_{i})P_{R}u(p_{j})$ and
$[ij]=\bar{u}(p_{i})P_{L}u(p_{j})$.

####  2.4.3 Matrix element including additional gluons

We can now use these results to form the analogue of equation (3) for the
process $gf_{b}\to H\cdots f_{b}$

$\displaystyle\begin{split}\overline{\left|\mathcal{M}_{{\tt HEJ}}^{gf_{b}\to
H\cdots\cdots
f_{b}}\right|}^{2}={}&\mathcal{B}_{H,f_{b}}(p_{a},p_{b},p_{1},p_{n},q_{j},q_{j+1})\\\
&\cdot\prod_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{n-2}\mathcal{V}(p_{a},p_{b},p_{a},p_{n},q_{i},q_{i+1})\\\
&\cdot\prod_{i=1}^{n-1}\mathcal{W}(q_{i},y_{i},y_{i+1}),\end{split}$ (18)

where the only differences to equation (3) are the Born-level function
$\mathcal{B}_{H,f_{b}}(p_{a},p_{b},p_{1},p_{n},q_{j},q_{j+1})$ given in
equation (12) and the third argument of the real-correction function
$\mathcal{V}$. We illustrate that this gives the correct behaviour in the MRK
limit in figure 3(c) for the processes $gu\to Hgu$, and in figure 3(d) we show
that we obtain the correct limiting behaviour for the NLL configuration $gu\to
Hug$.

##  3 Predictions and Comparison to Data

In this section we compare predictions for Higgs boson production in
association with one or more jets obtained with High Energy Jets to those of
fixed next-to-leading order perturbation theory and to experimental analyses.
The analyses are implemented in Rivet [29] and relate to data collected at the
LHC operated at both 13 TeV[30, 31] and 8 TeV [32].

###  3.1 Predictions

In our predictions, Sherpa [33] is used to generate leading-order events
through Comix [34] and Openloops [35] for $H+n$ jets, where $n=1,\dots,5$. We
include the exact dependence on the top-quark mass where available (i.e. for
$n=1,2$) and for higher multiplicities use the simpler results valid for an
infinite top mass. High-energy resummation is then applied using the method of
HEJ 2, described in detail in [23]. This takes the fixed-order events as input
and then adds all-order corrections (real and virtual) corresponding to each
Born phase space point. The resulting resummation events are reweighted by

$\displaystyle\frac{|{\cal M}_{\text{{\tt HEJ}}}(m_{t},m_{b})|^{2}}{|{\cal
M}_{\text{{\tt HEJ}, LO}}(m_{t},0)|^{2}}$ $\displaystyle\leq 2\text{ jets},$
(19) $\displaystyle\frac{|{\cal M}_{\text{{\tt
HEJ}}}(m_{t},m_{b})|^{2}}{|{\cal M}_{\text{{\tt HEJ}, LO}}(\infty,0)|^{2}}$
$\displaystyle>2\text{ jets}.$ (20)

${\cal M}_{\text{{\tt HEJ}}}(m_{t},m_{b})$ is the HEJ all-order matrix element
discussed in section 2, where we have indicated the dependence on the top-
quark mass $m_{t}$ and the bottom-quark mass $m_{b}$. ${\cal M}_{\text{{\tt
HEJ}, LO}}(m_{t},m_{b})$ denotes the leading-order truncation of the HEJ
matrix element. The $p_{T}$-sampling for the leading-order events used for the
matching extends slightly beyond the cuts used in the analysis, as required by
the mapping between the high-multiplicity $m$-body resummation phase space
point and the $n$-parton ($n<m$) phase space point of the matching. One way to
look at this is that the radiation produced by the resummation on top of the
fixed-order input modifies the momenta in the input, and the over-sampling is
needed in order for the full analysis-phase space of the resummation events to
be covered.

We also use Sherpa and Openloops to provide NLO 1-jet and 2-jet predictions in
the infinite top-quark mass limit without resummation, for comparisons with
HEJ and the experimental data. The cross sections presented from HEJ are
further matched to NLO by multiplying the predictions for the inclusive 1-jet
(or 2-jet) distributions by the ratio of the inclusive 1-jet (resp. 2-jet)
cross-section at NLO divided by the inclusive 1-jet (resp. 2-jet) cross-
section of HEJ expanded to NLO. This changes the normalisation of
distributions, and reduces the scale variation.

$\displaystyle\frac{d\sigma_{\text{HEJNLO}n\text{J}}}{d\mathcal{O}}={}$
$\displaystyle\frac{\sigma_{\text{NLO}n\text{J}}}{\sigma_{\text{{\tt
HEJ}}n\text{J}}}\frac{d\sigma_{\text{{\tt HEJ}}}}{d\mathcal{O}},$ (21)

where $\sigma_{\text{NLO}n\text{J}},n=1,2$ denotes the inclusive $n$-jet cross
section at NLO and $\sigma_{\text{{\tt HEJ}}n\text{J}}$ the HEJ prediction for
the inclusive $n$-jet cross section. Note that the components of the cross
section with exclusive three or more jets as predicted by HEJ are technically
matched only at Born level, but since they form part of the inclusive one or
two-jet observables, their contribution is scaled by the relevant ratio in eq.
(21).

We use the NNPDF30@NNLO [36] PDF set provided from the LHAPDF collaboration
[37] for HEJ and NLO predictions, with the central scale choice
$\mu_{F}=\mu_{R}=\max(m_{12},m_{H})$ (where $m_{12}$ is the invariant mass
between the two hardest jets, and set to $m_{12}=0$ for 1-jet events). In
order to gauge the scale dependence of the predictions the scales are varied
independently by a conventional factor of two, excluding combinations where
$\mu_{F}$ and $\mu_{R}$ differ by a factor of more than two. The coloured
regions in the figures below indicate the theoretical uncertainty envelope
formed by these scale variations.

We also investigated an alternative central scale choice
$\mu_{F}=\mu_{R}=H_{T}/2$. The predictions changed only minimally with this
scale compared to the custom scale choice above and so are not presented in
this study.

###  3.2 Predictions for 13 TeV and Comparison to Data

In this section we present predictions for a CMS analysis [30, 31] at a
centre-of-mass energy of $\sqrt{\hat{s}}=13$ TeV and for additional
distributions showcasing differences between HEJ and fixed order predictions
at NLO. The CMS study explored distributions for Higgs boson production (and
decay in the di-photon channel) both inclusively and in association with one
jet.

The baseline cuts related to the photons and the jets are listed in table 2
(see refs. [30, 31] for a full discussion). The pseudo-rapidity jet cuts are
specific to the observables studied and are listed in table 3. Jets are
reconstructed with the anti-$k_{T}$ [38] jet algorithm with $R=0.4$.

Description | Baseline cuts
---|---
Leading photon transverse momentum | $p_{T}(\gamma_{1})>30\text{ GeV}$
Subleading photon transverse momentum | $p_{T}(\gamma_{2})>18\text{ GeV}$
Diphoton invariant mass | $m_{\gamma\gamma}>90\text{ GeV}$
Pseudo-rapidity of the photons | $|\eta_{\gamma}|<2.5$
| excluding $1.4442<|\eta_{\gamma}|<1.566$
Ratio of harder photon $p_{T}$ to diphoton invariant mass | $p_{T}(\gamma_{1})/m_{\gamma\gamma}>\frac{1}{3}$
Ratio of softer photon $p_{T}$ to diphoton invariant mass | $p_{T}(\gamma_{2})/m_{\gamma\gamma}>\frac{1}{4}$
Photon isolation cut | $\text{Iso}^{\gamma}_{\text{gen}}<10\text{ GeV}$
Jet transverse momentum | $p_{T}(j)>30\text{ GeV}$
Table 2: Baseline photon and jet cuts of the 13 TeV analysis, following the CMS analysis of [30, 31]. $\text{Iso}^{\gamma}_{\text{gen}}$ denotes the sum of transverse energies of stable particles in a cone of radius $\Delta R$ = 0.3 around each photon. Observable | Pseudo-rapidity jet cut
---|---
Number of jets $N_{\text{jets}}$, figure 6(a) |  $|\eta_{j}|<2.5$ (all jets)
$|p_{T}^{j_{1}}|$, figure 6(b) |  $|\eta_{j_{1}}|<2.5$ (hardest jet) and $|\eta_{j}|<4.7$ (other jets)
$\min m_{ff}$, figure 7(a) |  $|\eta_{j_{1}}|<2.5$ (hardest jet) and $|\eta_{j}|<4.7$ (other jets)
$\max m_{ff}$, figure 7(b) |  $|\eta_{j_{1}}|<2.5$ (hardest jet) and $|\eta_{j}|<4.7$ (other jets)
Table 3: Pseudo-rapidity jet cuts used for the 13 TeV analysis observables
presented in this section, following the CMS analysis of [30, 31].

Figure 6(a) shows the exclusive number of jets where the 1-jet and 2-jet HEJ
predictions are rescaled as described in equation (21). The fixed-order
predictions are limited to 2 jets at NLO and 3 jets at LO, whereas HEJ allows
us to make predictions for the $\geq 4$-jet bin. The HEJ and NLO predictions
only describe $pp\to H+n$-jet processes via gluon fusion (GF) where the jets
consist of light quarks and gluons. The data includes a non-GF contribution
from electroweak VBF, $VH$ and $t\bar{t}H$ processes, labelled together as
$HX$ in the experimental papers. The size of this contribution is shown in the
analyses to be about 20% of the measured data values for the 2- and 3-jet bins
for example. Adding this to the HEJ and NLO results would increase both such
that the theoretical scale variation bands and the data uncertainty bands
overlap.

In figure 6(b), the transverse momentum of the first jet is shown. We have
compared to data from [30] here rather than [31] as it covers a larger range.
The non-GF contributions are not included in the NLO/HEJ predictions but are
included in the data, and are around 0.1 fb/GeV in the lower bins for example.
The discrepancy between NLO and HEJ predictions as the transverse momentum
increases is due to the resummation procedure, and has also been observed in
$W$+jets processes (see ref. [39]). The effect would be even more significant
for greater values of $p_{T}$, however the collected data does not probe this
region of phase-space. We have previously observed that a similar harder
$p_{T}$-spectrum seen in $H+\geq 2j$ processes in HEJ leads to a greater
sensitivity to the effects of using finite top and bottom quark masses [17].

The minimum rapidity separation between any two particles in the final state
is shown in figure 7(a). As the Higgs boson is one of these final states, this
is a 1-jet observable, so the NLO 1-jet predictions are shown for comparison
and the HEJ predictions are scaled by the ratio of the NLO to HEJ inclusive
1-jet rates. This observable is very sensitive to high-energy logarithmic
corrections, and as was observed in previous studies (see ref. [17]), the
effect of the resummation results in a significant lowering of the HEJ
prediction compared to fixed-order, by as much as 50% at large values. Figure
7(b) shows the maximum invariant mass between any two particles in the final
state. This is related to the high energy limit where all pairwise invariant
masses are taken to be large, but also includes situations where two or more
particles have a small invariant mass. The impact of the logarithmic
corrections is not as strong here, and the fixed-order and resummed
predictions agree within uncertainties.

(a)

(b)

Figure 6: (a) The exclusive number of jets compared to CMS data from [31], and
(b) the transverse momentum distribution of the leading jet compared to CMS
data [30]. Both analyses employ the cuts described in table 2.

(a)

(b)

Figure 7: High-energy sensitive inclusive 1-jet distributions: (a) the minimum
rapidity separation between any two outgoing particles (Higgs boson or jets)
and (b) the maximum invariant mass between any two outgoing particles (Higgs
or jets). HEJ results are rescaled by the inclusive cross section ratio
$\sigma_{\text{NLO1J}}/\sigma_{\text{{\tt HEJ} 1J}}$.

###  3.3 Predictions for 8 TeV and Comparison to Data

We now present predictions for an ATLAS analysis [32] at a centre-of-mass
energy of $\sqrt{\hat{s}}=8$ TeV as implemented in Rivet[29]. We list the
relevant experimental cuts used in this analysis in table 4, the complete list
being available in the experimental publication. As in the experimental
analysis, the jets are reconstructed with the anti-$k_{T}$ algorithm with a
radius parameter of $R=0.4$. This study explored the inclusive and
differential cross-sections for Higgs boson production in the diphoton decay
channel. For our purposes, we select the observables which correspond to Higgs
boson production plus at least one jet, where our predictions are applicable.

Description | Baseline cuts
---|---
Photon transverse momentum | $p_{T}(\gamma)>25\text{ GeV}$
Diphoton invariant mass | $105\text{ GeV}<m_{\gamma\gamma}<160\text{ GeV}$
Pseudo-rapidity of the photons |  $|\eta_{\gamma}|<2.37$ excluding $1.37<|\eta_{\gamma}|<1.56$
Ratio of harder photon $p_{T}$ to diphoton invariant mass | $p_{T}(\gamma_{1})/m_{\gamma\gamma}>0.35$
Ratio of softer photon $p_{T}$ to diphoton invariant mass | $p_{T}(\gamma_{2})/m_{\gamma\gamma}>0.25$
Photon isolation cut | $\text{Iso}^{\gamma}_{\text{gen}}<14\text{ GeV}$
Jet transverse momentum | $p_{T}(j)>30\text{ GeV}$
Jet rapidity | $|y_{j}|<4.4$
Table 4: Baseline cuts of the 8 TeV analysis, following the ATLAS analysis of
[32]. $\text{Iso}^{\gamma}_{\text{gen}}$ denotes the sum of transverse
energies of stable particles in a cone of radius $\Delta R$ = 0.4 around each
photon.

We divide our results into 1-jet observables, i.e. containing at least one
jet, where the new components of HEJ as detailed in section 2.3 can be tested,
and 2-jet observables.

####  3.3.1 $H+\geq 1j$

In figure 8(a), we show the exclusive number of jets. As was evidenced at 13
TeV, the differences between fixed-order and resummed predictions are limited
after the inclusive cross sections are rescaled. As in the previous
comparison, the experimental data points here include a non-GF contribution of
a few fb. Adding this to the NLO and HEJ predictions would increase
predictions for the 1- and 2-jet rates such that the bands for the theoretical
scale variance and data uncertainty bands overlap. The predictions in the
$\geq 3$-jet bin remain slightly below data even with the extra contribution.

In figure 8(b), the rapidity of the leading jet is displayed: the discrepancy
between the fixed-order and the resummed predictions increases as the rapidity
of jet attains large values. This is a High-Energy effect as opposed to a
finite quark mass effect. Indeed, the corrections in $\hat{s}/t$ are
particularly sizeable in this region of phase-space, and previous studies
(see. ref [17]) showed little dependence on the inclusion of the finite quark
mass effects on this observable. However, this is not the case for the
transverse momentum of the Higgs boson of figure 8(c). The finite quark mass
effects and the resummation lead to a hardening of the high-$p_{T}$ tail of
the Higgs boson, which would be even more dramatic had that region been
probed. Due to the probed phase-space region of $p_{T}<140\,\text{GeV}$, HEJ
and fixed-order predictions for the hardest jet transverse momentum of figure
8(d) remain close together and difficult to disentangle.

The High-Energy sensitive observables of figure 9 behave in a similar fashion
to those at 13 TeV (figure 7) for the reasons explained in section 3.2. Here,
there is nearly a factor of two difference between NLO and HEJ at large values
of min $\Delta y_{ff}$.

(a) Number of jets

(b) Leading jet rapidity

(c) $p_{T}$ of the Higgs boson, $N_{\text{jets}}=1$

(d) Leading jet $p_{T}$, $N_{\text{jets}}=1$

Figure 8: (8(a)): Number of jets (exclusive). (8(b)): Leading jet rapidity.
(8(c)): Higgs boson transverse momentum in the 1-jet bin. (8(d)): Leading jet
transverse momentum in the 1-jet bin. The 1-jet HEJ predictions are rescaled
by the inclusive cross section ratio $\sigma_{\text{NLO1J}}/\sigma_{\text{{\tt
HEJ} 1J}}$ while the HEJ predictions of the 2 and 3-jet bins of (8(a)) are
rescaled by $\sigma_{\text{NLO2J}}/\sigma_{\text{{\tt HEJ} 2J}}$.

(a) $\text{min }\Delta y_{ff}$

(b) $\text{max }m_{ff}$

Figure 9: High-energy sensitive 1-jet distributions. (9(a)): minimum rapidity
separation between any two outgoing particles (Higgs boson or jets). (9(b)):
maximum invariant mass between any two outgoing particles (Higgs boson or
jets). HEJ results are rescaled by the inclusive cross section ratio
$\sigma_{\text{NLO1J}}/\sigma_{\text{{\tt HEJ} 1J}}$.

####  3.3.2 $H+\geq 2j$

We now turn to a range of 2-jet observables, displayed in figures 10 and 11.
As in the previous sections, the data contains a non-GF component absent from
NLO and HEJ predictions. Adding this increases both sets of predictions closer
to data, in some cases such that the scale variance bands overlap with the
uncertainty bands on the data.

Globally, the impact of the resummation on High-Energy sensitive observables
is to lower the predictions from fixed-order approaches, as can be seen in
large dijet rapidity separation in figure 10(a), large rapidity values of the
second hardest jet in figure 10(b) and at large dijet invariant mass in figure
10(c). As expected, the resummation procedure has little impact on the
observables dependent on the azimuthal degrees of freedom: the azimuthal angle
difference between the leading two jets of figure 10(d) and the azimuthal
angle difference between the diphoton and the leading dijets systems depicted
in figure 11(a), expect perhaps at values close to $\pi$ (that is when the
systems are back-to-back).

As previously observed, the combination of the inclusion of corrections in
$\hat{s}/t$ and the finite quark mass effects tend to harden the tail of the
transverse momenta distributions compared to fixed order predictions. This is
apparent in the description of the third-leading jet transverse momentum of
figure 11(b), but also in the transverse momentum of the diphoton-dijet system
of figure 11(c). Although the Higgs transverse momentum seems to be
independent of the effect of the resummation, it is conjectured that values of
$p^{\gamma\gamma}_{\perp}$ above 200 GeV would lead to a disparity between the
two approaches.

(a) Dijet rapidity separation

(b) Subleading jet rapidity

(c) Invariant dijet mass

(d) $\Delta\phi$ between the leading 2 jets

Figure 10: (10(a)): Dijet rapidity separation. (10(b)): Subleading jet
rapidity. (10(c)): dijet invariant mass. (10(d)): Azimuthal angle difference
between the leading 2 jets. All 2-jet HEJ predictions are rescaled by the
inclusive cross section ratio $\sigma_{\text{NLO2J}}/\sigma_{\text{{\tt HEJ}
2J}}$.

(a) $\Delta\phi$ between dijet and diphoton systems

(b) Third-leading jet $p_{T}$

(c) Higgs $p_{T}$ with $N_{\text{jets}}=2$

(d) $p_{T}$ of diphoton-dijet system

Figure 11: (11(a)): Azimuthal angle difference between dijet and diphoton
objects. (11(b)): Transverse momentum of the third-leading jet. (11(c)):
Transverse momentum of the Higgs boson in the $\geq$ 2-jet bin. (11(d)):
Transverse momentum of the Higgs plus dijet object :
$(p_{H}+p_{j_{1}}+p_{j_{2}})_{\perp}$. All 2 and 3-jet HEJ predictions are
rescaled by the inclusive cross section ratio
$\sigma_{\text{NLO2J}}/\sigma_{\text{{\tt HEJ} 2J}}$.

##  4 Conclusions

In this paper we have presented an alternative description of $pp\to H+\geq
1j$, which is accurate to leading logarithms in $\hat{s}/p_{T}^{2}$ (LL). We
have outlined the structure of a LL-accurate amplitude in the HEJ formalism,
and described the calculation of the necessary new components in section 2.
One big advantage of the approach is that it maintains full dependence on the
finite top and bottom quark masses in the couplings of the Higgs boson to
gluons for _any_ number of jets, which quickly exceeds the multiplicities
currently calculated at even leading order. The new pieces allow LL
resummation in $\hat{s}/p_{T}^{2}$ to an inclusive 1-jet process for the first
time in the HEJ framework.

We have then compared the resummed predictions to fixed-order predictions and
to LHC data in section 3, and discussed the impact of the logarithmic
corrections. We find the impact of the resummation is seen at large jet
transverse momenta. The resummed results give a harder $p_{T}$-spectrum
compared to NLO, which in turn leads to a greater dependence on finite quark
masses in the coupling. We also observe a large suppression compared to NLO at
large values of rapidity separation between all pairs of final state particles
(i.e. between any two of the Higgs boson and jets). This can be as much of a
factor of two and lies significantly outwith the uncertainty bands on the two
predictions. Other observables, e.g. azimuthal angles, are less sensitive to
these logarithmic corrections.

Looking forward to analyses of LHC Run 3 data, our results suggest that the
inclusion of finite quark masses for higher jet multiplicities _and_ of
logarithmic corrections in $\hat{s}/p_{T}^{2}$ will be important in the
comparison to data.

## Acknowledgements

We are grateful to the other members of the HEJ collaboration for useful and
helpful discussions throughout this work. We are pleased to acknowledge
funding from the UK Science and Technology Facilities Council (under grant
number ST/T506047/1 for HH), the Royal Society and the ERC Starting Grant
715049 “QCDforfuture”. The predictions presented in section 3 were produced
using resources from PhenoGrid which is part of the GridPP Collaboration [40,
41]. AP acknowledges support by the National Science Foundation under Grant
No. PHY 2210161. For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

##  5 Appendix - NLO reweighting factors

In table 5, we give the value of the NLO reweighting factors as described in
equation (21) for both the 8 TeV [32] and 13 TeV [30, 31] analyses.

Analysis | 8 TeV | 13 TeV
---|---|---
Scale | $\mu_{F},\mu_{R}$ | $(\mu_{F},\mu_{R})/2$ | $2(\mu_{F},\mu_{R})$ | $\mu_{F},\mu_{R}$ | $(\mu_{F},\mu_{R})/2$ | $2(\mu_{F},\mu_{R})$
1J factor | 1.87 | 1.54 | 2.15 | 1.59 | 1.30 | 1.84
2J factor | 1.98 | 1.48 | 2.40 | 1.62 | 1.19 | 2.00
Table 5: NLO Reweighting factors with $\mu_{F}=\mu_{R}=\max(m_{12},m_{H})$.

In the 8 TeV analysis, the inclusive 1-jet and 2-jet cross-sections are
calculated from the rapidity of the hardest and second hardest jet histograms
respectively, figures 8(b) and 10(b).

In the 13 TeV analysis, the inclusive 1-jet and 2-jet cross-sections are
obtained from the appropriate $N_{\text{jets}}$ bins, figure 6(a). Note that
for at least 2 jets, this plot requires central jets only, but it is valid to
use it as it is the only plot we present for the 2-jet observables. If more
inclusive cross-sections are considered, say from the rapidity of the second
hardest jet histogram over all the experimental range, then the 2-jet
reweighting factor would be further away from the 1-jet value (1.79 for the
central scale instead of 1.62).

## References

  * [1] M. Czakon, R. V. Harlander, J. Klappert and M. Niggetiedt, _Exact Top-Quark Mass Dependence in Hadronic Higgs Production_ , _Phys. Rev. Lett._ 127 (2021) 162002, [2105.04436].
  * [2] S. P. Jones, M. Kerner and G. Luisoni, _Next-to-Leading-Order QCD Corrections to Higgs Boson Plus Jet Production with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 120 (2018) 162001, [1802.00349].
  * [3] X. Chen, A. Huss, S. P. Jones, M. Kerner, J. N. Lang, J. M. Lindert et al., _Top-quark mass effects in H+jet and H+2 jets production_ , _JHEP_ 03 (2022) 096, [2110.06953].
  * [4] V. Del Duca, W. Kilgore, C. Oleari, C. Schmidt and D. Zeppenfeld, _Higgs \+ 2 jets via gluon fusion_ , _Phys. Rev. Lett._ 87 (2001) 122001, [hep-ph/0105129].
  * [5] V. Del Duca, W. Kilgore, C. Oleari, C. Schmidt and D. Zeppenfeld, _Gluon fusion contributions to H + 2 jet production_ , _Nucl. Phys. B_ 616 (2001) 367–399, [hep-ph/0108030].
  * [6] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, _Higgs Boson Gluon-Fusion Production in QCD at Three Loops_ , _Phys. Rev. Lett._ 114 (2015) 212001, [1503.06056].
  * [7] F. Dulat, B. Mistlberger and A. Pelloni, _Differential Higgs production at N 3LO beyond threshold_, _JHEP_ 01 (2018) 145, [1710.03016].
  * [8] B. Mistlberger, _Higgs boson production at hadron colliders at N 3LO in QCD_, _JHEP_ 05 (2018) 028, [1802.00833].
  * [9] L. Cieri, X. Chen, T. Gehrmann, E. W. N. Glover and A. Huss, _Higgs boson production at the LHC using the $q_{T}$ subtraction formalism at N3LO QCD_, _JHEP_ 02 (2019) 096, [1807.11501].
  * [10] X. Chen, T. Gehrmann, E. W. N. Glover, A. Huss, B. Mistlberger and A. Pelloni, _Fully Differential Higgs Boson Production to Third Order in QCD_ , _Phys. Rev. Lett._ 127 (2021) 072002, [2102.07607].
  * [11] R. Boughezal, F. Caola, K. Melnikov, F. Petriello and M. Schulze, _Higgs boson production in association with a jet at next-to-next-to-leading order in perturbative QCD_ , _JHEP_ 06 (2013) 072, [1302.6216].
  * [12] X. Chen, T. Gehrmann, E. W. N. Glover and M. Jaquier, _Precise QCD predictions for the production of Higgs + jet final states_ , _Phys. Lett. B_ 740 (2015) 147–150, [1408.5325].
  * [13] R. Boughezal, F. Caola, K. Melnikov, F. Petriello and M. Schulze, _Higgs boson production in association with a jet at next-to-next-to-leading order_ , _Phys. Rev. Lett._ 115 (2015) 082003, [1504.07922].
  * [14] J. M. Campbell, R. K. Ellis and G. Zanderighi, _Next-to-Leading order Higgs + 2 jet production via gluon fusion_ , _JHEP_ 10 (2006) 028, [hep-ph/0608194].
  * [15] J. M. Campbell, R. K. Ellis and C. Williams, _Hadronic Production of a Higgs Boson and Two Jets at Next-to-Leading Order_ , _Phys. Rev. D_ 81 (2010) 074023, [1001.4495].
  * [16] V. Del Duca, W. Kilgore, C. Oleari, C. R. Schmidt and D. Zeppenfeld, _Kinematical limits on Higgs boson production via gluon fusion in association with jets_ , _Phys. Rev. D_ 67 (2003) 073003, [hep-ph/0301013].
  * [17] J. R. Andersen, J. D. Cockburn, M. Heil, A. Maier and J. M. Smillie, _Finite Quark-Mass Effects in Higgs Boson Production with Dijets at Large Energies_ , _JHEP_ 04 (2019) 127, [1812.08072].
  * [18] B.-W. Xiao and F. Yuan, _BFKL and Sudakov Resummation in Higgs Boson Plus Jet Production with Large Rapidity Separation_ , _Phys. Lett. B_ 782 (2018) 28–33, [1801.05478].
  * [19] F. G. Celiberto, D. Y. Ivanov, M. M. A. Mohammed and A. Papa, _High-energy resummed distributions for the inclusive Higgs-plus-jet production at the LHC_ , _Eur. Phys. J. C_ 81 (2021) 293, [2008.00501].
  * [20] J. R. Andersen and J. M. Smillie, _Constructing All-Order Corrections to Multi-Jet Rates_ , _JHEP_ 1001 (2010) 039, [0908.2786].
  * [21] J. R. Andersen and J. M. Smillie, _The Factorisation of the t-channel Pole in Quark-Gluon Scattering_ , _Phys.Rev._ D81 (2010) 114021, [0910.5113].
  * [22] J. R. Andersen and J. M. Smillie, _Multiple Jets at the LHC with High Energy Jets_ , _JHEP_ 1106 (2011) 010, [1101.5394].
  * [23] J. R. Andersen, T. Hapola, M. Heil, A. Maier and J. M. Smillie, _Higgs-boson plus Dijets: Higher-Order Matching for High-Energy Predictions_ , _JHEP_ 08 (2018) 090, [1805.04446].
  * [24] V. S. Fadin, R. Fiore, M. G. Kozlov and A. V. Reznichenko, _Proof of the multi-Regge form of QCD amplitudes with gluon exchanges in the NLA_ , _Phys. Lett. B_ 639 (2006) 74–81, [hep-ph/0602006].
  * [25] J. R. Andersen, T. Hapola, A. Maier and J. M. Smillie, _Higgs Boson Plus Dijets: Higher Order Corrections_ , _JHEP_ 09 (2017) 065, [1706.01002].
  * [26] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * [27] S. Dawson and R. P. Kauffman, _Higgs boson plus multi - jet rates at the SSC_ , _Phys. Rev. Lett._ 68 (1992) 2273–2276.
  * [28] R. P. Kauffman, S. V. Desai and D. Risal, _Production of a Higgs boson plus two jets in hadronic collisions_ , _Phys. Rev. D_ 55 (1997) 4005–4015, [hep-ph/9610541].
  * [29] C. Bierlich et al., _Robust Independent Validation of Experiment and Theory: Rivet version 3_ , _SciPost Phys._ 8 (2020) 026, [1912.05451].
  * [30] CMS collaboration, A. M. Sirunyan et al., _Measurement of inclusive and differential Higgs boson production cross sections in the diphoton decay channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 01 (2019) 183, [1807.03825].
  * [31] CMS collaboration, _Measurement of the Higgs boson inclusive and differential fiducial production cross sections in the diphoton decay channel with pp collisions at $\sqrt{s}$ = 13 TeV_, 2208.12279.
  * [32] ATLAS collaboration, G. Aad et al., _Measurements of fiducial and differential cross sections for Higgs boson production in the diphoton decay channel at $\sqrt{s}=8$ TeV with ATLAS_, _JHEP_ 09 (2014) 112, [1407.4222].
  * [33] Sherpa collaboration, E. Bothmann et al., _Event Generation with Sherpa 2.2_ , _SciPost Phys._ 7 (2019) 034, [1905.09127].
  * [34] T. Gleisberg and S. Höche, _Comix, a new matrix element generator_ , _Journal of High Energy Physics_ 2008 (Dec, 2008) 039–039.
  * [35] F. Buccioni, S. Pozzorini and M. Zoller, _On-the-fly reduction of open loops_ , _Eur. Phys. J. C_ 78 (2018) 70, [1710.11452].
  * [36] R. D. Ball, , V. Bertone, S. Carrazza, C. S. Deans, L. D. Debbio et al., _Parton distributions for the LHC run II_ , _Journal of High Energy Physics_ 2015 (apr, 2015) .
  * [37] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _The European Physical Journal C_ 75 (mar, 2015) .
  * [38] M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063, [0802.1189].
  * [39] J. R. Andersen, J. A. Black, H. M. Brooks, E. P. Byrne, A. Maier and J. M. Smillie, _Combined subleading high-energy logarithms and NLO accuracy for W production in association with multiple jets_ , _JHEP_ 04 (2021) 105, [2012.10310].
  * [40] The GridPP Collaboration, _GridPP: Development of the UK Computing Grid for Particle Physics_ , _J. Phys. G_ 32 (2006) N1–N20.
  * [41] D. Britton et al., _GridPP: the UK grid for particle physics_ , _Phil. Trans. R. Soc. A_ 367 (2009) 2447–2457.


j}\widetilde{R}_{j}^{(h)^{-1}}\widetilde{H}^{(h)}_{i\to
j}\right\\}_{h=1}^{k}+K_{i}\boldsymbol{F}_{i}K_{i}^{\top}\right)K_{i},$ (19)

where $\oplus$ is the direct sum operator,
$\boldsymbol{F}_{i}=\oplus\\{\boldsymbol{A}_{i}(\boldsymbol{\ell})\\}_{\boldsymbol{\ell}\in{\cal
S}_{i}}$, and after letting
$x_{j}(\boldsymbol{\ell})=\boldsymbol{\lambda}_{[j,:]}\boldsymbol{v}(\boldsymbol{\ell})$,
we compute
$\boldsymbol{A}_{i}(\boldsymbol{\ell})=-\sum_{j=1}^{q}\boldsymbol{\lambda}_{[j,:]}^{\top}\boldsymbol{\lambda}_{[j,:]}E\left[\frac{\delta^{2}}{\delta^{2}x_{j}(\boldsymbol{\ell})}\log
dF_{j}(y_{j}(\boldsymbol{\ell})\,|\,\boldsymbol{v}(\boldsymbol{\ell}),\boldsymbol{\lambda}_{[j,:]},\beta_{j},\gamma_{j})\right]$.

##  Appendix E Applications supplement

In all our applications, all methods we implement are allowed to use up to 16
CPU threads in a workstation with 128GB memory and an AMD Ryzen 9 5950X CPU. R
package meshed (v.0.2) allows to set the number of OpenMP (Dagum and Menon, ,
1998) threads, whereas Hmsc takes advantage of parallelization via BLAS when
performing expensive operations (e.g., chol($\cdot$) ). The R-INLA package
used to implement SPDE-INLA methods can similarly take advantage of
multithreaded operations.

###  E.1 Model outputs in the bivariate count data application

We report model outputs from the SPDE-INLA and QMGP-MALA models for the
bivariate count data application of Section 5.1.

Figure 16: Output from fitting a SPDE model via INLA. Top row: recovered log-
intensity and predictions for both outcomes. Bottom row: width of posterior
credible intervals about log-intensity, and residual log-intensity. Figure 17:
Output from fitting a coregionalized QMGP via MALA. Top row: recovered log-
intensity and predictions for both outcomes. Bottom row: width of posterior
credible intervals about log-intensity, and residual log-intensity.

###  E.2 Simulating 750 datasets with bivariate count outcomes

We provide additional details on the 750 simulated datasets of Figure 8 in the
main article. We generate Poisson data on a $50\times 50$ regular grid, for a
total of 2500 observations for $y_{j}(\boldsymbol{\ell})\sim
Pois(\exp\\{\eta_{j}(\boldsymbol{\ell})\\})$ where
$\boldsymbol{\eta}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{v}(\boldsymbol{\ell})$
and $\boldsymbol{v}(\cdot)$ is a bivariate GP with independent Matérn
correlations with $\nu_{j}=1/2$ for $j=1,2$ and $\phi_{2}=2.5$. We choose
$\phi_{1}\in\\{2.5,12.5,25\\}$. We introduce missing values at $1/5$ of
spatial locations independently for each outcome. We fix the $2\times 2$
loading matrix via

$\boldsymbol{\Lambda}=\begin{bmatrix}\lambda_{11}&\lambda_{12}\\\
\lambda_{21}&\lambda_{22}\end{bmatrix}=\texttt{chol}\left(\begin{bmatrix}\lambda_{1}&0\\\
0&1\end{bmatrix}\cdot\begin{bmatrix}1&\rho\\\
\rho&1\end{bmatrix}\cdot\begin{bmatrix}\lambda_{1}&0\\\
0&1\end{bmatrix}\right),$

which implies $\lambda_{11}=\lambda_{1}$, $\lambda_{12}=0$, and $\lambda_{21}$
and $\lambda_{22}$ are such the latent correlation between the first and
second margin is exactly $\rho$. We choose
$\lambda_{1}\in\\{\frac{\sqrt{2}}{2},2\\}$ and
$\rho\in\\{-.9,-0.25,0,0.65,.9\\}$. We generate 25 datasets for all
combination of values of $\phi_{1}$, $\lambda_{1}$ and $\rho$. All QMGP models
implemented in Section 5.1 use the GriPS method for sampling covariance
parameters (Peruzzi et al., , 2021). Since the data are gridded, GriPS amounts
to a parameter expansion strategy similar to the one used in Ghosh and Dunson,
(2009). We report comparative results of the same analyses with/without GriPS
in Figure 18. GriPS does not affect the relative performance of the tested
methodologies for sampling QMGPs; however, it improves sampler efficiency and
uncertainty quantification about the linear predictor. We note here that GriPS
is less effective in higher dimensional multivariate cases and we did not
observe advantages in the community ecology applications of Section 6 in the
main article.

Figure 18: Estimation of $\rho$, its ESS/s, and coverage of intervals about
$\eta_{j,\text{test}}$, with GriPS. Methods without GriPS are denoted as
“-def”.

###  E.3 Latent process sampler efficiency in multi-type data

We use the same setup of Section E.2, choosing $\lambda_{1}\in\\{2.5,12.5\\}$,
for each of the following pairs of outcome types:
$\\{(\text{Gaussian},\text{Poisson}),(\text{Neg.
Binomial},\text{Binomial}),(\text{Neg. Binomial},\text{Poisson})\\}$, for a
total of 1500 datasets, of which 500 include a Binomial or Gaussian outcome,
1000 include a Poisson or Neg. Binomial outcome. We fix all unknowns
($\boldsymbol{\Lambda}$, $\phi_{1}$, $\phi_{2}$) to their true values except
for the latent process, as we are interested in purely comparing sampling
methods for $\boldsymbol{w}(\cdot)$. For each dataset, we calculate ESS/s for
samples of $\boldsymbol{w}(\boldsymbol{\ell}_{i})$,
$\boldsymbol{\ell}_{i}=1,\dots n$. After computing the median ESS/s as a
summary efficiency measure for each dataset, we compute the mean of this
measure over all datasets. Efficiency summary results are reported in Table 5.
Comparisons based on RMSPE and coverage about
$\eta_{j}(\boldsymbol{\ell}),j=1,2$ are in Table 6

Method | Binomial | Gaussian | Negative Binomial | Poisson
---|---|---|---|---
MALA | 1.15 | 8.35 | 1.31 | 2.81
NUTS | 0.25 | 2.15 | 0.32 | 0.68
SiMPA | 4.26 | 17.65 | 9.03 | 8.90
SM-MALA | 1.81 | 7.85 | 3.45 | 3.59
Table 5: Efficiency in posterior sampling of $\boldsymbol{w}(\cdot)$, in terms
of ESS/s, for different types of outcome in the bivariate synthetic data
application with multi-type outcomes.

Method | RMSPE | Covg. 95%
---|---|---
Binomial | Gaussian | Neg. Bin. | Poisson | Binomial | Gaussian | Neg. Bin. | Poisson
MALA | 0.450 | 0.332 | 14.206 | 1.198 | 0.897 | 0.942 | 0.899 | 0.930
NUTS | 0.449 | 0.328 | 13.932 | 1.195 | 0.932 | 0.944 | 0.934 | 0.941
SiMPA | 0.449 | 0.327 | 13.923 | 1.194 | 0.944 | 0.948 | 0.947 | 0.947
SM-MALA | 0.449 | 0.327 | 13.971 | 1.212 | 0.943 | 0.948 | 0.939 | 0.940

Table 6: RMSPE in predicting the outcomes
$\boldsymbol{y}_{j}(\boldsymbol{\ell})$ on the test set, and empirical
coverage of 95% credible intervals about $\eta_{j}(\boldsymbol{\ell})$,
$j=1,2$, for different types of outcome in the bivariate synthetic data
application with multi-type outcomes.

## References

  * Andrieu and Thoms, (2008) Andrieu, C. and Thoms, J. (2008).  A tutorial on adaptive MCMC.  Statistics and Computing, 18:343–373.  doi:10.1007/s11222-008-9110-y.
  * Apanasovich and Genton, (2010) Apanasovich, T. V. and Genton, M. G. (2010).  Cross-covariance functions for multivariate random fields based on latent dimensions.  Biometrika, 97:15–30.  doi:10.1093/biomet/asp078.
  * Atchadé, (2006) Atchadé, Y. F. (2006).  An adaptive version for the Metropolis adjusted Langevin algorithm with a truncated drift.  Methodology and Computing in Applied Probability, 8:235–254.  doi:10.1007/s11009-006-8550-0.
  * Banerjee, (2017) Banerjee, S. (2017).  High-dimensional Bayesian geostatistics.  Bayesian Analysis, 12(2):583–614.  doi:10.1214/17-BA1056R.
  * Banerjee, (2020) Banerjee, S. (2020).  Modeling massive spatial datasets using a conjugate Bayesian linear modeling framework.  Spatial Statistics, 37:100417.  doi:10.1016/j.spasta.2020.100417.
  * Banerjee et al., (2010) Banerjee, S., Finley, A. O., Waldmann, P., and Ericsson, T. (2010).  Hierarchical spatial process models for multiple traits in large genetic trials.  Journal of American Statistical Association, 105(490):506–521.  doi:10.1198/jasa.2009.ap09068.
  * Banerjee et al., (2008) Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008).  Gaussian predictive process models for large spatial data sets.  Journal of the Royal Statistical Society, Series B, 70:825–848.  doi:10.1111/j.1467-9868.2008.00663.x.
  * Betancourt, (2018) Betancourt, M. (2018).  A conceptual introduction to Hamiltonian Monte Carlo.  arXiv:1701.02434.
  * Bhattacharya and Dunson, (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse Bayesian infinite factor models.  Biometrika, 98(2):291–306.  doi:10.1093/biomet/asr013.
  * Blomstedt et al., (2019) Blomstedt, P., Parente Paiva Mesquita, D., Lintusaari, J., Sivula, T., Corander, J., and Kaski, S. (2019).  Meta-analysis of Bayesian analyses.  arXiv:1904.04484.
  * Bradley et al., (2018) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2018).  Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion).  Bayesian Analysis, 13(1):253–310.  doi:10.1214/17-BA1069.
  * Bradley et al., (2019) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2019).  Bayesian hierarchical models with conjugate full-conditional distributions for dependent data from the natural exponential family.  Journal of the American Statistical Association.  doi:10.1080/01621459.2019.1677471.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, Marcus Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1).  doi:10.18637/jss.v076.i01.
  * Chen et al., (2020) Chen, Z., Wang, B., and Gorban, A. N. (2020).  Multivariate Gaussian and Student-t process regression for multi-output prediction.  Neural Computing and Applications, 32:3005–3028.  doi:10.1007/s00521-019-04687-8.
  * Cressie and Johannesson, (2008) Cressie, N. and Johannesson, G. (2008).  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society, Series B, 70:209–226.  doi:10.1111/j.1467-9868.2007.00633.x.
  * Dagum and Menon, (1998) Dagum, L. and Menon, R. (1998).  OpenMP: an industry standard api for shared-memory programming.  Computational Science & Engineering, IEEE, 5(1):46–55.
  * (17) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016a).  Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association, 111:800–812.  doi:10.1080/01621459.2015.1044091.
  * (18) Datta, A., Banerjee, S., Finley, A. O., Hamm, N. A. S., and Schaap, M. (2016b).  Nonseparable dynamic nearest neighbor Gaussian process models for large spatio-temporal data with an application to particulate matter analysis.  The Annals of Applied Statistics, 10:1286–1316.  doi:10.1214/16-AOAS931.
  * Dey et al., (2021) Dey, D., Datta, A., and Banerjee, S. (2021).  Graphical Gaussian process models for highly multivariate spatial data.  Biometrika.  in press. doi:doi.org/10.1093/biomet/asab061.
  * Duane et al., (1987) Duane, S., A.D., K., Pendleton, B. J., and Roweth, D. (1987).  Hybrid Monte Carlo.  Physics Letters B, 195:216–222.
  * Dunson and Johndrow, (2020) Dunson, D. and Johndrow, J. E. (2020).  The Hastings algorithm at fifty.  Biometrika, 107(1):1–23.  doi:10.1093/biomet/asz066.
  * Finley et al., (2008) Finley, A. O., Banerjee, S., Ek, A. R., and McRoberts, R. E. (2008).  Bayesian multivariate process modeling for prediction of forest attributes.  Journal of Agricultural, Biological, and Environmental Statistics, 13:60.  doi:10.1198/108571108X273160.
  * Finley et al., (2019) Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019).  Efficient algorithms for Bayesian nearest neighbor Gaussian processes.  Journal of Computational and Graphical Statistics, 28:401–414.  doi:10.1080/10618600.2018.1537924.
  * Furrer et al., (2006) Furrer, R., Genton, M. G., and Nychka, D. (2006).  Covariance tapering for interpolation of large spatial datasets.  Journal of Computational and Graphical Statistics, 15:502–523.  doi:10.1198/106186006X132178.
  * Genton and Kleiber, (2015) Genton, M. G. and Kleiber, W. (2015).  Cross-covariance functions for multivariate geostatistics.  Statistical Science, 30:147–163.  doi:10.1214/14-STS487.
  * Ghosh and Dunson, (2009) Ghosh, J. and Dunson, D. B. (2009).  Default prior distributions and efficient posterior computation in Bayesian factor analysis.  Journal of Computational and Graphical Statistics, 18(2):306–320.  doi:10.1198/jcgs.2009.07145.
  * Girolami and Calderhead, (2011) Girolami, M. and Calderhead, B. (2011).  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B, 73(2):123–214.  doi:10.1111/j.1467-9868.2010.00765.x.
  * Gramacy and Apley, (2015) Gramacy, R. B. and Apley, D. W. (2015).  Local Gaussian process approximation for large computer experiments.  Journal of Computational and Graphical Statistics, 24(2):561–578.  doi:10.1080/10618600.2014.914442.
  * Guhaniyogi and Banerjee, (2018) Guhaniyogi, R. and Banerjee, S. (2018).  Meta-kriging: Scalable Bayesian modeling and inference for massive spatial datasets.  Technometrics, 60(4):430–444.  doi:10.1080/00401706.2018.1437474.
  * Haario et al., (2001) Haario, H., Saksman, E., and Tamminen, J. (2001).  An adaptive Metropolis algorithm.  Bernoulli, 7(2):223–242.  doi:10.2307/3318737.
  * Heaton et al., (2019) Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit-Mangion, A. (2019).  A case study competition among methods for analyzing large spatial data.  Journal of Agricultural, Biological and Environmental Statistics, 24(3):398–425.  doi:10.1007/s13253-018-00348-w.
  * Hoffman and Gelman, (2014) Hoffman, M. D. and Gelman, A. (2014).  The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15:351–1381.
  * Jin et al., (2021) Jin, B., Peruzzi, M., and Dunson, D. B. (2021).  Bag of DAGs: Flexible & scalable modeling of spatiotemporal dependence.  arXiv:2112.11870.
  * Johndrow et al., (2020) Johndrow, J. E., Pillai, N. S., and Smith, A. (2020).  No free lunch for approximate mcmc.  arXiv:2010.125147.
  * Jönsson et al., (2010) Jönsson, A. M., Eklundh, L., Hellström, M., Bärring, L., and Jönsson, P. (2010).  Annual changes in MODIS vegetation indices of Swedish coniferous forests in relation to snow dynamics and tree phenology.  Remote Sensing of Environment, 114:2719–2730.  doi:10.1016/j.rse.2010.06.005.
  * Jurek and Katzfuss, (2020) Jurek, M. and Katzfuss, M. (2020).  Hierarchical sparse Cholesky decomposition with applications to high-dimensional spatio-temporal filtering.  arXiv:2006.16901.
  * Katzfuss, (2017) Katzfuss, M. (2017).  A multi-resolution approximation for massive spatial datasets.  Journal of the American Statistical Association, 112:201–214.  doi:10.1080/01621459.2015.1123632.
  * Katzfuss and Guinness, (2021) Katzfuss, M. and Guinness, J. (2021).  A general framework for Vecchia approximations of Gaussian processes.  Statistical Science, 36(1):124–141.  doi:10.1214/19-STS755.
  * Kaufman et al., (2008) Kaufman, C. G., Schervish, M. J., and Nychka, D. W. (2008).  Covariance tapering for likelihood-based estimation in large spatial data sets.  Journal of the American Statistical Association, 103:1545–1555.  doi:10.1198/016214508000000959.
  * Lindgren et al., (2011) Lindgren, F., Rue, H., and Lindström, J. (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B, 73:423–498.  doi:10.1111/j.1467-9868.2011.00777.x.
  * Marshall and Roberts, (2012) Marshall, T. and Roberts, G. (2012).  An adaptive approach to Langevin MCMC.  Statistics and Computing, 22:1041–1057.  doi:10.1007/s11222-011-9276-6.
  * Matheron, (1982) Matheron, G. (1982).  Pour une analyse krigeante des données régionalisées.  Technical report N.732, Centre de Géostatistique.
  * Mesquita et al., (2020) Mesquita, D., Blomstedt, P., and Kaski, S. (2020).  Embarrassingly parallel MCMC using deep invertible transformations.  In Adams, R. P. and Gogate, V., editors, Proceedings of Machine Learning Research, volume 115, pages 1244–1252, Tel Aviv, Israel. PMLR.  http://proceedings.mlr.press/v115/mesquita20a.html.
  * Neal, (2011) Neal, R. M. (2011).  MCMC using Hamiltonian dynamics.  In Brooks, S., Gelman, A., Jones, G. L., and Meng, X.-L., editors, Handbook of Markov Chain Monte Carlo. CRC Press, New York.  doi:10.1201/b10905.
  * Neiswanger et al., (2014) Neiswanger, W., Wang, C., and Xing, E. P. (2014).  Asymptotically exact, embarrassingly parallel MCMC.  In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI’14, page 623–632, Arlington, Virginia, USA. AUAI Press.
  * Nemeth and Sherlock, (2018) Nemeth, C. and Sherlock, C. (2018).  Merging mcmc subposteriors through gaussian-process approximations.  Bayesian Analysis, 13(2):507–530.  doi:doi = "10.1214/17-BA1063".
  * Peruzzi et al., (2021) Peruzzi, M., Banerjee, S., Dunson, D. B., and Finley, A. O. (2021).  Grid-Parametrize-Split (GriPS) for improved scalable inference in spatial big data analysis.  arXiv:2101.03579.
  * Peruzzi et al., (2020) Peruzzi, M., Banerjee, S., and Finley, A. O. (2020).  Highly scalable Bayesian geostatistical modeling via meshed Gaussian processes on partitioned domains.  Journal of the American Statistical Association.  in press. doi:10.1080/01621459.2020.1833889.
  * Peruzzi and Dunson, (2021) Peruzzi, M. and Dunson, D. B. (2021).  Spatial multivariate trees for big data Bayesian regression.  Journal of Machine Learning Research, accepted. arXiv:2012.00943.
  * Roberts and Rosenthal, (2007) Roberts, G. O. and Rosenthal, J. S. (2007).  Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms.  Journal of Applied Probability, 44:458–475.  doi:10.1239/jap/1183667414.
  * Roberts and Stramer, (2002) Roberts, G. O. and Stramer, O. (2002).  Langevin diffusions and Metropolis-Hastings algorithms.  Methodology And Computing In Applied Probability, 4:337–357.  doi:10.1023/A:1023562417138.
  * Roberts and Tweedie, (1996) Roberts, G. O. and Tweedie, R. L. (1996).  Exponential convergence of Langevin distributions and their discrete approximations.  Bernoulli, 2(4):341–363.
  * Rue and Held, (2005) Rue, H. and Held, L. (2005).  Gaussian Markov Random Fields: Theory and Applications.  Chapman & Hall/CRC.  doi:10.1007/978-3-642-20192-9.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested laplace approximations.  Journal of the Royal Statistical Society: Series B, 71:319–392.  doi:10.1111/j.1467-9868.2008.00700.x.
  * Sang and Huang, (2012) Sang, H. and Huang, J. Z. (2012).  A full scale approximation of covariance functions for large spatial data sets.  Journal of the Royal Statistical Society, Series B, 74:111–132.  doi:10.1111/j.1467-9868.2011.01007.x.
  * Schmidt and Gelfand, (2003) Schmidt, A. M. and Gelfand, A. E. (2003).  A Bayesian coregionalization approach for multivariate pollutant data.  Journal of Geophysical Research, 108:D24.  doi:10.1029/2002JD002905.
  * Sengupta and Cressie, (2013) Sengupta, A. and Cressie, N. (2013).  Hierarchical statistical modeling of big spatial datasets using the exponential family of distributions.  Spatial Statistics.  doi:10.1016/j.spasta.2013.02.002.
  * Shah et al., (2014) Shah, A., Wilson, A. G., and Ghahramani, Z. (2014).  Student-t processes as alternatives to Gaussian processes.  In Proceedings of the 17th International Conference on Artificial Intelligence and Statistics (AISTATS).
  * Shirota et al., (2019) Shirota, S., Finley, A. O., Cook, B. D., and Banerjee, S. (2019).  Conjugate nearest neighbor Gaussian process models for efficient statistical interpolation of large spatial data.  arXiv:1907.10109.
  * Stein, (2014) Stein, M. L. (2014).  Limitations on low rank approximations for covariance matrices of spatial data.  Spatial Statistics, 8:1–19.  doi:doi:10.1016/j.spasta.2013.06.003.
  * Stein et al., (2004) Stein, M. L., Chi, Z., and Welty, L. J. (2004).  Approximating likelihoods for large spatial data sets.  Journal of the Royal Statistical Society, Series B, 66:275–296.  doi:10.1046/j.1369-7412.2003.05512.x.
  * Sun et al., (2011) Sun, Y., Li, B., and Genton, M. (2011).  Geostatistics for large datasets.  In Montero, J., Porcu, E., and Schlather, M., editors, Advances and Challenges in Space-time Modelling of Natural Events, pages 55–77. Springer-Verlag, Berlin Heidelberg.  doi:10.1007/978-3-642-17086-7.
  * Taylor and Diggle, (2014) Taylor, B. M. and Diggle, P. J. (2014).  INLA or MCMC? a tutorial and comparative evaluation for spatial prediction in log-Gaussian Cox processes.  Journal of Statistical Computation and Simulation, 84(10):2266–2284.  doi:10.1080/00949655.2013.788653.
  * Tikhonov et al., (2020) Tikhonov, G., Opedal, O. H., Abrego, N., Lehikoinen, A., de Jonge, M. M. J., Oksanen, J., and Ovaskainen, O. (2020).  Joint species distribution modelling with the R-package Hmsc.  Methods in Ecology and Evolution, 11(3):442–447.  doi:10.1111/2041-210X.13345.
  * Vecchia, (1988) Vecchia, A. V. (1988).  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society, Series B, 50:297–312.  doi:10.1111/j.2517-6161.1988.tb01729.x.
  * Vihola, (2012) Vihola, M. (2012).  Robust adaptive Metropolis algorithm with coerced acceptance rate.  Statistics and Computing, 22:997–1008.  doi:10.1007/s11222-011-9269-5.
  * Wackernagel, (2003) Wackernagel, H. (2003).  Multivariate Geostatistics: An Introduction with Applications.  Springer, Berlin.  doi:10.1007/978-3-662-05294-5.
  * Walker et al., (1993) Walker, D. A., Halfpenny, J. C., Walker, M. D., and Wessman, C. A. (1993).  Long-term studies of snow-vegetation interactions.  BioScience, 43(5):287–301.  doi:10.2307/1312061.
  * (69) Wang, K., Zhang, L., Qiu, Y., Ji, L., Tian, F., Wang, C., and Wang, Z. (2015a).  Snow effects on alpine vegetation in the Qinghai-Tibetan plateau.  International Journal of Digital Earth, 8(1):58–75.  doi:10.1080/17538947.2013.848946.
  * Wang and Dunson, (2014) Wang, X. and Dunson, D. B. (2014).  Parallelizing MCMC via Weierstrass sampler.  arXiv:1312.4605.
  * (71) Wang, X., Guo, F., Heller, K. A., and Dunson, D. B. (2015b).  Parallelizing MCMC with random partition trees.  In Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume 1, NIPS’15, page 451–459, Cambridge, MA, USA. MIT Press.  arXiv:1506.03164.
  * Xie et al., (2020) Xie, J., Jonas, T., Rixen, C., de Jong, R., Garonna, I., Notarnicola, C., Asam, S., Schaepman, M. E., and Kneubühler, M. (2020).  Land surface phenology and greenness in Alpine grasslands driven by seasonal snow and meteorological factors.  Science of The Total Environment, 725:138380.  doi:10.1016/j.scitotenv.2020.138380.
  * Zanella and Roberts, (2021) Zanella, G. and Roberts, G. (2021).  Multilevel linear models, gibbs samplers and multigrid decompositions.  Bayesian Analysis.  doi:10.1214/20-BA1242.
  * Zhang and Banerjee, (2021) Zhang, L. and Banerjee, S. (2021).  Spatial factor modeling: A Bayesian matrix-normal approach for misaligned data.  Biometrics.  doi:10.1111/biom.13452.
  * Zilber and Katzfuss, (2020) Zilber, D. and Katzfuss, M. (2020).  Vecchia-Laplace approximations of generalized Gaussian processes for big non-Gaussian spatial data.  arXiv:1906.07828.


#  nrCascadeSim \- A simulation tool for nuclear recoil cascades resulting
from neutron capture

A.N. Villano  Kitty Harris  Staci Brown

††margin:  DOI: https://doi.org/DOI unavailablehttps://doi.orgDOI unavailable
Software • N/Ahttps://doi.orgReview • NO_REPOSITORYhttps://doi.orgRepository •
DOI unavailablehttps://doi.orgArchive Editor:
https://example.comhttps://doi.orgPending Editor
Reviewers: • https://github.com/Pending Reviewershttps://doi.org@Pending
Reviewers Submitted: N/A
Published: N/A License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

Neutron capture-induced nuclear recoils have emerged as an important tool for
detector calibrations in direct dark matter detection and coherent elastic
neutrino-nucleus scattering (CE$\mathrm{\nu}$NS).

nrCascadeSim is a C++ command-line tool for generating simulation data for
energy deposits resulting from neutron capture on pure materials. Presently,
capture events within silicon, germanium, neon, and argon are supported. While
the software was developed for solid state detector calibration, it can be
used for any application which requires simulated neutron capture-induced
nuclear recoil data.

A “cascade” occurs when a neutron becomes part of a nucleus. The neutron can
be captured to one of many discrete energy levels, or states; if the energy
level is nonzero (not the ground state), then the state will eventually change
so that it is zero. This can happen either all at once or in multiple steps —
that is, the captured neutron may go from its state to the ground state, or it
may go to another state with lower energy that is not the ground state
(provided that one exists). The cascade refers to the particular “path” of
energy levels that a captured neutron takes to get to the ground state from
the neutron separation energy. Currently, the code assumes that the neutrons
that enter the nuclear system have zero kinetic energy; this is a good
approximation for thermal neutrons because 0.0254 eV (the average kinetic
energy of a thermal neutron) is small compared to most nuclear recoil energy
scales, making it negligible.

nrCascadeSim models many of these cascades at once and saves the energies
along with other useful data to a single file. The output file is a ROOT file
(Brun & Rademakers, 1997).

## Models Used

When modeling deposits from neutron capture events, we want to look at the
recoil of the atom as a result of these cascades. To determine how much energy
is deposited, we must track how much the atom slows down between steps of the
cascade as well as how each nuclear state change affects the atom’s kinetic
energy. nrCascadeSim assumes a constant deceleration that results from the
atom colliding with other nearby electrons and nuclei. This means that it must
simulate, along with the steps of the cascade, the time between each state —
to calculate how much the atom slows down. And it must also simulate the angle
between the atom’s momentum before a decay and the momentum boost (gamma ray)
resulting from the decay — to calculate the resulting momenta. The time
between steps is simulated as an exponential random variable based on the
state’s half-life, and the angle is simulated as having an isotropic
distribution. Cascade selection is weighted by isotope abundance (Rosman &
Taylor, 1998; Sonzogni & Shu, 2020) and cross-section as well as the
probability of the energy level. In existing levelfiles, energy levels are
derived from (Islam et al., 1991) for germanium and from (Raman et al., 1992)
for silicon.

The above process models the recoil energies, and the output gives both the
total recoil energy for a cascade as well as the energy per step. For some
applications, this may be the desired output, or the user may already have a
particular process they will use for converting this energy to what they wish
to measure. However, we also include, for convenience, the ionization yield
and ionization energy of these recoils. Ionization yield is a fraction that,
when multiplied by the energy, gives the ionization energy, and ionization
energy is the amount of energy that would be read out if an otherwise
equivalent electron recoil were to occur. This calculation is useful because
many solid-state detectors read out the ionization energy for nuclear recoils.
This ionization yield assumes the Lindhard model (Lindhard et al., 1963).

Figure 1 compares the normalized frequencies of ionization energies from the
Lindhard (Lindhard et al., 1963) model with the Sorensen (Sorensen, 2015)
yield model, which is applied after the simulation using Python, and applies
detector resolution models to both. This figure demonstrates one example of
user-applied analysis utilizing the energy deposits at each step instead of
the ionization energy.

Figure 1: An overlaid histogram showing an example use case in which points
are generated and then multiple yield models and resolutions are applied. The
“Small Res (1/5)” histograms have Gaussians with 1/5 of the width of their
counterparts.

## Statement of Need

The goal of this software is to simplify the computation of the nuclear recoil
spectrum following neutron capture for a variety of applications. These
include nuclear recoil calibrations for dark matter direct detection and
coherent neutrino detection (CE$\mathrm{\nu}$NS). In these cases as the
particle detection has become more sensitive (detectors having a lower energy
threshold) it is now possible to use the capture-induced nuclear recoil events
for detector calibrations. Additionally, thermalized neutrons will provide
large backgrounds that have heretofore not been modeled. The key roadblock to
studying these scenarios is the complexity of calculating the nuclear recoil
spectrum.

nrCascadeSim addresses this need by allowing users to generate nuclear recoil
simulations that reflect a variety of single-element detector setups. The
energy levels that the recoiling nuclei may pass between and their respective
lifetimes are customizable, and multiple isotopes of the same element can be
present within the same simulation. Pre-defined energy level files exist for
silicon and germanium, which take into account the natural abundance data of
each isotope in (Rosman & Taylor, 1998) and (Sonzogni & Shu, 2020). Output
values include energy deposits at each step along each individual cascade,
total kinetic energy deposits, and ionization energy deposits.

## State of the Field

While there are tools, such as the open-source GEANT4 (Agostinelli et al.,
2003) framework, that allow users to simulate neutron capture, existing tools
are not built specifically for neutron capture-based nuclear recoils as
nrCascadeSim is and therefore use some underlying assumptions that
nrCascadeSim does not. The main approximation often used in GEANT4 that we
avoid in nrCascadeSim is that all recoils decay directly to the ground state.
While this works for some applications, it is necessary to be more precise
when an accurate spectrum of neutron capture-based recoils is needed for
analyses such as calibration or background subtraction. Figure 2 shows a
comparison for the energy deposits produced by Geant4 for natural silicon
compared with those produced by nrCascadeSim. The figure does not include any
instrumentation resolution and shows a highly prominent peak around 1.25 keV
recoil energy (coming from capture on 29Si directly to the ground state)
whereas the nrCascadeSim shows another direct-to-ground contribution (from
capture on 28Si) at around 1.0 keV recoil energy and generally far more
“spread out” recoils coming from two- or more step cascades.

Figure 2: An overlaid histogram showing how the Geant4 v10.7.3 energy deposits
compare with those from nrCascadeSim for natural silicon.

Recently, the power of the neutron capture-induced events has been
acknowledged in the CE$\mathrm{\nu}$NS field (Thulliez & others, 2021). That
initial study, however, used the FIFRELIN code (Litaize & Serot, 2010), which
was originally developed for modeling fission fragments and has been updated
to use statistical models of gamma emission for the purpose of modeling
fission-fragment deexcitation (Litaize et al., 2014). nrCascadeSim takes the
complementary approach of beginning with small to medium-sized nuclei and
modeling the cascades in more exact detail. The goal is for the code to be
extended to heavier nuclei but still using this detailed approach.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics (HEP) under Award Number DE-
SC0021364.

## References

Agostinelli, S., Allison, J., Amako, K., Apostolakis, J., Araujo, H., Arce,
P., Asai, M., Axen, D., Banerjee, S., Barrand, G., Behner, F., Bellagamba, L.,
Boudreau, J., Broglia, L., Brunengo, A., Burkhardt, H., Chauvie, S., Chuma,
J., Chytracek, R., & Cooperman, G. (2003). GEANT4 - a simulation toolkit.
_Nuclear Instruments and Methods in Physics Research A_ , _506_ , 250–303.
https://doi.org/10.1016/S0168-9002(03)01368-8

Brun, R., & Rademakers, F. (1997). ROOT - an object oriented data analysis
framework. _Nucl. Inst. & Meth. In Phys. Res. A_, _389_ , 81–86. See also
"ROOT" [software], Release v6.22/00, 02/07/2020, (No DOI available for this
version, closest version linked). https://doi.org/10.5281/zenodo.3895852

Islam, M. A., Kennett, T. J., & Prestwich, W. V. (1991). Radiative strength
functions of germanium from thermal neutron capture. _Physical Review C_ ,
_43_ , 1086–1098. https://doi.org/10.1103/physrevc.43.1086

Lindhard, J., Nielsen, V., & Scharff, M. (1963). Integral equations governing
radiation effects. _Mat. Fys. Medd. Dan. Vid. Selsk._ , _33_ , 1–41.

Litaize, O., & Serot, O. (2010). Investigation of phenomenological models for
the monte carlo simulation of the prompt fission neutron and $\gamma$
emission. _Phys. Rev. C_ , _82_ , 054616\.
https://doi.org/10.1103/PhysRevC.82.054616

Litaize, O., Serot, O., & Berge, L. (2014). https://www.oecd-
nea.org/science/meetings/pnd22/presentations/1-LITAIZE.pdf

Raman, S., Jurney, E. T., & Lynn, J. E. (1992). Thermal-neutron capture by
silicon isotopes. _Physical Review C_ , _46_ , 972–983.
https://doi.org/10.1103/PhysRevC.46.972

Rosman, K. J. R., & Taylor, P. D. P. (1998). Isotopic compositions of the
elements 1997. In _Pure & Appl. Chem._ (Vol. 70, pp. 217–235). International
Union of Pure; Applied Chemistry. https://doi.org/10.1515/iupac.70.0026

Sonzogni, A., & Shu, B. (2020). _Nudat 2_. Brookhaven National Laboratory.
https://www.nndc.bnl.gov/nudat2/

Sorensen, P. (2015). Atomic limits in the search for galactic dark matter.
_Phys. Rev. D_ , _91_. https://doi.org/10.1103/PhysRevD.91.083509

Thulliez, L., & others. (2021). _Calibration of nuclear recoils at the 100 eV
scale using neutron capture_. _16_(07), P07032.
https://doi.org/10.1088/1748-0221/16/07/p07032