CERN-TH-2023-229

disorder: Deep inelastic scattering at high orders

Alexander Karlberg

CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland

## Abstract

We present a Fortran 77/95 code capable of computing QCD corrections in deep
inelastic scattering (DIS). The code uses the Projection-to-Born method to
combine an exclusive DIS + 1 jet calculation with the inclusive DIS structure
functions, thereby obtaining fully differential DIS predictions at
$\mathcal{O}(\alpha_{\mathrm{s}}^{2})$. The code is lightweight and fast, and
yet includes the most common functionalities found in typical perturbative QCD
programs, like automatic renormalisation and factorisation scale
uncertainties, options to run and combine multiple seeds, and interfaces to
fastjet and LHAPDF. Due to the underlying exclusive DIS + 1 jet code, the
program also provides stable results in the infrared, relevant for extracting
logarithmic coefficients for analytic resummations. As a by-product the code
provides access to the DIS structure functions and (reduced) cross sections up
to $\mathcal{O}(\alpha_{\mathrm{s}}^{3})$.

https://github.com/alexanderkarlberg/disorder

###### Contents

  1. 1 Introduction
  2. 2 Basics of the DIS process
    1. 2.1 A note on the phase space
    2. 2.2 Applying P2B
  3. 3 Running disorder
    1. 3.1 Compiling and prerequisites
      1. 3.1.1 The analysis framework
    2. 3.2 Inclusive mode
    3. 3.3 P2B mode
    4. 3.4 Validating the code
  4. 4 Benchmarks and results
    1. 4.1 Inclusive results
    2. 4.2 Exclusive laboratory frame results
  5. 5 Conclusion

##  1 Introduction

Deep inelastic scattering (DIS) is arguably one of the best understood
processes in perturbative QCD. It is one of few processes for which there
exists an exact factorisation theorem [1, 2], and the massless unpolarised
hard perturbative coefficients have been computed through an impressive three
loops [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] with progress on the four
loop results presented recently in Ref. [16]. Together with the three-loop
results for the splitting functions [17, 18, 19] and the four-loop
$\beta$-function [20, 21] for the running coupling this allows for the
determination of the proton structure functions at the next-to-next-to-next-
to-leading order (N3LO).111Technically the four-loop splitting functions are
needed to claim this accuracy. Recent progress in determining those and the
three-loop massive operator matrix elements can be found in Refs. [22, 23, 24,
25, 26, 27, 28] and [29, 30, 31, 32, 33, 34, 35] respectively.
Phenomenologically the first few moments which have been computed therein,
with some exact pieces also being known, are however in all likelihood enough
to claim N3LO.  The structure functions can be combined with an exclusive
next-to-next-to-leading order (NNLO) DIS + 1 jet computation to obtain fully
differential N3LO predictions in DIS, as was done by the NNLOjet collaboration
[36, 37]. This computation uses the Projection-to-Born (P2B) method, which was
first introduced in the context of NNLO Vector Boson Fusion (VBF) production
[38].

Despite this impressive theoretical progress, there are few publically
available computer codes from which one can obtain fast and reliable high-
order differential cross section predictions for DIS. This paper and
associated Fortran code seeks to address that.

One of the advantages of the P2B method is that it is rather agnostic towards
the details of the underlying exclusive computation, and it is hence possible
to “upgrade” existing fixed-order codes to one order higher. Historically the
most successful fixed-order codes have been MEPJET [39], disent [40],
disaster++ [41], and nlojet++ [42] which are all next-to-leading order (NLO)
accurate. Discrepancies between disent and fixed-order coefficients from
analytical resummation were initially observed in Refs.[43, 44], and were only
recently understood to be due to a bug in one of the dipole terms for the
gluon channel in disent [45, 46].

With the bug fixed we can use disent as the underlying exclusive NLO code
together with the NNLO structure functions from hoppet [47, 48].222The
combination of disent and hoppet naturally leads to a Fortran code. For a C++
alternative, based on publically available codes, one could have started from
either disaster++ or nlojet++ and used the structure functions as implemented
in APFEL++ [49, 50] Advantages of using disent are its well-known efficiency
and numerical stability. The resulting program is dubbed disorder keeping in
the spirit of previous names for fixed-order DIS codes.333dispatch [44]
deserves an honorable mentioning in this context. Since the program relies on
hoppet for the structure functions, the inclusive DIS cross section can
actually be obtained one order higher, as the structure functions were
implemented at N3LO already in the context of the proVBFH codes [38, 51, 52,
53] and are publically available since v1.3.0 of hoppet [54]444hoppet
implements all known approximate and exact ingredients at this order - see
footnote 1.. The disorder program therefore has two more or less separate use
cases: 1. The computation of fully exclusive NNLO photon-mediated massless DIS
2. The computation of inclusive neutral- (NC) and charged-current (CC)
massless DIS at N3LO accuracy.

The program is designed to be user-friendly with a simple command line
interface. It prints cross sections and all important run-parameters to the
screen and disk, allowing a user to acquire cross sections with very little
effort. It comes with an interface to fastjet [55] and LHAPDF [56] and uses
the histogramming package from the POWHEG-BOX [57] for easy analysis. The code
can also compute renormalisation and factorisation scale uncertainties on-the-
fly. For inclusive cross sections the code provides results in a matter of
seconds555The majority of this time is taken up by the structure function
initialisation inside hoppet rather than the integration of the cross section.
even at N3LO with integration uncertainties typically below the permille
level. For exclusive quantities the code can be run on a laptop at NLO and
depending on the analysis, and the laptop, even at NNLO.

It should be pointed out that the code provides little theoretical advance on
its own. As described above, the structure functions themselves have been
known for a little while, and the P2B method has already been applied to DIS,
even one order higher than here.666It has also been applied in polarised DIS
in Refs. [45, 46, 58, 59, 60]. It is however the author’s opinion that with
the renewed interest in DIS due to the upcoming Electron-Ion-Collider (EIC)
[61], well-maintained public code is extremely valuable for both the
experimental and theoretical communities, and that providing documentation in
the form of this article will enable the wide use of the code.

The paper is structured as follows: In section 2 we review the DIS process and
kinematics as implemented in disorder. In section 3 we provide details on how
to run disorder and in section 4 we show a few results from the program.
Finally we conclude in section 5.

##  2 Basics of the DIS process

In this section we first give some standard definitions for kinematics and
cross sections in DIS. Along the way we specify the conventions that are used
in the disorder code, and provide some details on the P2B method as applied
here.

At leading order (LO) the DIS process is the scattering of a massless
(anti-)quark $q$ off a massless (anti-)lepton $l$ via the exchange of a photon
or electroweak gauge boson $V$ of virtuality $Q^{2}$. Denoting the external
four-momenta by $k_{i}$ (incoming lepton), $k_{f}$ (outgoing lepton), $p_{i}$
(incoming quark), and $p_{f}$ (outgoing quark) we can define the Lorentz
invariant DIS variables $x$, $Q^{2}$, and $y$, given by

$\displaystyle Q^{2}=-q^{2}=-(k_{i}-k_{f})^{2},\qquad x=\frac{Q^{2}}{2P\cdot
q},\qquad y=\frac{P\cdot q}{P\cdot k_{i}}=\frac{p_{i}\cdot q}{p_{i}\cdot
k_{i}},$ (1)

where $P$ is the hadron four-momentum. As can be seen these kinematics are
fully specified by the hadron and lepton momenta. This is also true beyond LO.
DIS is most often analysed in the Breit-frame which is specified777For the
explicit transformation between lab- and Breit-frames we follow Appendix 7.11
in Ref. [62]. by requiring that $2x\vec{P}+\vec{q}=0$. In this frame the
mediated vector boson has zero energy component and is anti-aligned with the
incoming parton. Explicitly in disorder the Breit-frame at LO is

$\displaystyle k_{i}$
$\displaystyle=\frac{Q}{2}\left(\frac{2-y}{y},\frac{2\sqrt{1-y}}{y},0,-1\right),\quad
p_{i}=\frac{Q}{2}\left(1,0,0,+1\right)$ (2) $\displaystyle k_{f}$
$\displaystyle=\frac{Q}{2}\left(\frac{2-y}{y},\frac{2\sqrt{1-y}}{y},0,+1\right),\quad
p_{f}=\frac{Q}{2}\left(1,0,0,-1\right)$

where four-momenta are given as $(E,p_{x},p_{y},p_{z})$. The resulting vector
$q$ is hence given by $(0,0,0,-Q)$. In the laboratory frame we align the
parton with the positive $z$-axis and hence the lepton with the negative

$\displaystyle\tilde{k}_{i}$
$\displaystyle=E_{\mathrm{l}}\left(1,0,0,-1\right),\quad\tilde{p}_{i}=xE_{\mathrm{h}}\left(1,0,0,+1\right)$
(3)

where $E_{\mathrm{h}}$ is the energy of the incoming hadron. Using the
definitions of eq. (1) one finds that in the lab frame

$\tilde{q}=\left(y(E_{\mathrm{l}}-xE_{\mathrm{h}}),-Q\sqrt{1-y},0,-y(E_{\mathrm{l}}-xE_{\mathrm{h}})\right)$
(4)

and the outgoing momenta then simply follow from momentum conservation

$\tilde{k}_{f}=\tilde{k}_{i}-\tilde{q},\quad\tilde{p}_{f}=\tilde{p}_{i}+\tilde{q}\,.$
(5)

The inclusive cross section for DIS can be split into a NC contribution, from
$e^{\pm}p\to e^{\pm}+X$ scattering, and a CC contribution from
$e^{\pm}p\to\nu+X$888Here $\nu=\\{\nu_{e},\bar{\nu}_{e}\\}$ for $e^{-}$ and
$e^{+}$ respectively. One can of course also consider incoming neutrinos which
does not change the discussion here. scattering. The unpolarised NC cross
section can be written as

$\frac{\mathrm{d}\sigma_{\mathrm{NC}}^{\pm}}{\mathrm{d}x\mathrm{d}Q^{2}}=\frac{2\pi\alpha^{2}}{xQ^{4}}\left[y_{+}F_{2}^{\mathrm{NC}}\mp
y_{-}xF_{3}^{\mathrm{NC}}-y^{2}F_{L}^{\mathrm{NC}}\right],\,$ (6)

where $y_{\pm}=1\pm(1-y)^{2}$, $\alpha$ is the fine structure constant and
$F_{i}^{\mathrm{NC}}$ can be expressed in terms of the usual proton structure
functions

$\displaystyle F_{i}^{\mathrm{NC}}$
$\displaystyle=F_{i}^{\gamma}-v_{e}\Gamma_{\gamma Z}F_{i}^{\gamma
Z}+(v_{e}^{2}+a_{e}^{2})\Gamma_{Z}F_{i}^{Z},\,\quad i=2,L$ (7) $\displaystyle
F_{3}^{\mathrm{NC}}$ $\displaystyle=-a_{e}\Gamma_{\gamma Z}F_{3}^{\gamma
Z}+2v_{e}a_{e}\Gamma_{Z}F_{3}^{Z}\,.$

Here $v_{e}=-\frac{1}{2}+2\sin^{2}\theta_{W}$ is the vector and
$a_{e}=\frac{1}{2}$ is the axial-vector coupling, $M_{Z}$ is the Z boson mass,
$\Gamma_{\gamma Z}=\frac{Q^{2}}{\sin^{2}2\theta_{W}(Q^{2}+M_{Z}^{2})}$, and
$\Gamma_{Z}=\Gamma_{\gamma Z}^{2}$. $\theta_{W}$ is the weak mixing angle. In
disorder the electroweak parameters are fixed by $\alpha$, $M_{W}$, and
$M_{Z}$ through the tree-level relations

$\displaystyle\sin^{2}\theta_{W}=1-\frac{M_{W}^{2}}{M_{Z}^{2}},\qquad
G_{F}=\frac{\pi\alpha}{\sqrt{2}m_{W}^{2}\sin^{2}\theta_{W}},$ (8)

with $G_{F}$ the Fermi constant. Similarly we define the unpolarised CC cross
section as

$\frac{\mathrm{d}\sigma_{\mathrm{CC}}^{\pm}}{\mathrm{d}x\mathrm{d}Q^{2}}=\frac{\pi\alpha^{2}}{8\sin^{4}\theta_{W}x}\left[\frac{1}{M_{W}^{2}+Q^{2}}\right]^{2}\left[y_{+}F_{2}^{\mathrm{CC}}\mp
y_{-}xF_{3}^{\mathrm{CC}}-y^{2}F_{L}^{\mathrm{CC}}\right],\,$ (9)

where $M_{W}$ is the mass of $W$ boson and the CC structure functions are now
simply given by the $W$ ones

$F_{2}^{\mathrm{CC}}=F_{2}^{W^{\pm}},\quad
F_{L}^{\mathrm{CC}}=F_{L}^{W^{\pm}},\quad
F_{3}^{\mathrm{CC}}=F_{3}^{W^{\pm}}\,.$ (10)

The exact definitions of all proton structure functions inside hoppet up to
N3LO can be found in Refs. [47, 54, 48]. The tabulation of the structure
functions inside hoppet, as used in disorder, has a relative numerical
precision of $10^{-4}$ for most values of $x$, which in turn limits the
accuracy which can be obtained to the same order.

In DIS it is also customary to define the dimensionless _reduced_ NC and CC
cross sections by [63]

$\displaystyle\tilde{\sigma}_{\mathrm{NC}}^{\pm}(x,Q^{2})$
$\displaystyle=\frac{xQ^{4}}{2\pi\alpha^{2}}\frac{1}{y_{+}}\frac{\mathrm{d}\sigma_{\mathrm{NC}}^{\pm}}{\mathrm{d}x\mathrm{d}Q^{2}}$
(11) $\displaystyle\tilde{\sigma}_{\mathrm{CC}}^{\pm}(x,Q^{2})$
$\displaystyle=\frac{8\sin^{4}\theta_{W}x}{\pi\alpha^{2}}\left[M_{W}^{2}+Q^{2}\right]^{2}\frac{\mathrm{d}\sigma_{\mathrm{CC}}^{\pm}}{\mathrm{d}x\mathrm{d}Q^{2}}\,.$

disorder provides direct access to all the cross sections in eqs. (6), (9),
and (11) at N3LO accuracy. In principle one can also access the NC and CC
structure functions, although they are currently only computed as an
intermediate step to construct the cross sections and are not accessible in
the user analysis.

###  2.1 A note on the phase space

Since the Born phase space is fully constrained by specifying any two of the
three DIS variables $x$, $y$, and $Q^{2}$ there exists more than one double-
differential cross section. The convention in disorder is to always return
$\frac{\mathrm{d}\sigma^{2}}{\mathrm{d}x\mathrm{d}Q^{2}}$ (in
pb/$\,\mathrm{GeV}^{2}$) regardless of which of the two variables are fixed.
The user can easily convert to
$\frac{\mathrm{d}\sigma^{2}}{\mathrm{d}x\mathrm{d}y}$ or
$\frac{\mathrm{d}\sigma^{2}}{\mathrm{d}y\mathrm{d}Q^{2}}$ by supplying a
factor $\frac{Q^{2}}{y}$ or $\frac{x}{y}$ respectively. If only one of $x$,
$y$, or $Q^{2}$ is specified then the appropriate single differential cross
section is returned, i.e. $\frac{\mathrm{d}\sigma}{\mathrm{d}x}$ in pb,
$\frac{\mathrm{d}\sigma}{\mathrm{d}y}$ in pb, or
$\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}}$ in pb/$\,\mathrm{GeV}^{2}$,
integrated over the other two variables. If none of the variables are fixed
then the total cross section integrated over $x$, $y$, and $Q^{2}$ is returned
in pb.

###  2.2 Applying P2B

The structure functions are by definition inclusive in all radiation and can
therefore only provide predictions for quantities, like the inclusive cross
sections, which depend on the Born kinematics of eq. (1) only. If we instead
evaluate the structure functions on an observable sensitive to emissions, e.g.
the transverse momentum of the hardest jet in the laboratory frame, we see
that this will not give the right answer, as the real emissions are not
included with their correct kinematics.

The P2B method lifts this restriction by effectively replacing the Born-
kinematics real-emission contributions in the structure functions with the
correct kinematics ones. In practice whenever disent returns an event with
some weight, we bin it once according to the true kinematics, and again
projecting the kinematics to the underlying Born changing the sign of the
weight. This last term, upon integration, will exactly cancel the real
contribution in the structure functions, whereas the first term will provide
the correct real matrix element.

The projections themselves are trivial because, as outlined above, the Born
kinematics are fully specified at all orders by the lepton (and proton)
momenta, as given in eqs. (2)–(5). Hence computing the projections adds very
little computational effort to the cross section calculation.

At this point it is worth reminding the reader that disent only includes the
photon-mediated NC, and hence disorder only provides exclusive predictions for
this channel. It should in principle be possible to extend the code to include
Z-mediation (including interferences) and CC, but we leave this for future
work.

In principle P2B could be applied to the disent code without any major
modifications besides the bug fix mentioned in the introduction. However a few
significant modifications were introduced to allow for a more flexible
integration of disent into disorder. The version of disent that we include in
disorder is based on the version which can be found in dispatch [44], and
which already included some minor modifications. In addition to the
modifications present there, we have also introduced $\alpha$ as an input
parameter, whereas before it was fixed to $1/137$. A number of parameters can
also now be set on the command line as described below. Besides coding the
actual interface to disent, the biggest modification introduced in disorder is
that the KPFUNS subroutine now returns an array of weights corresponding to
varying the factorisation scale by a factor two up and down. This allows for
much faster evaluation of scale uncertainties compared to running the program
three separate times. Additionally disent can now also use any scale defined
in terms of $x$, $y$, and $Q$, and not just some multiple of $Q$. We provide a
few different scale choices as documented below, but more can easily be
implemented.

##  3 Running disorder

In this section we give instructions on how to compile and run disorder,
giving a few examples of the use of the most common command line arguments.
The code itself can be obtained from

https://github.com/alexanderkarlberg/disorder

###  3.1 Compiling and prerequisites

A user should start by inspecting README.md. To compile disorder both hoppet
(v1.3.0 or later) and LHAPDF (tested with v6.5.4) [56] have to be installed on
the machine. If both are installed in a location in the $PATH it is enough to
run

⬇

mkdir build && cd build

cmake ..

make [-j]

from the main directory. This will create an executable disorder along with
two auxiliary executables, mergedata and getpdfuncert. For non-standard
installation of hoppet and LHAPDF the paths can be specified like this

⬇

cmake -DHOPPET_CONFIG=/path/to/hoppet-config -DLHAPDF_CONFIG=/path/to/lhapdf-
config

where the path should include the config-file itself (i.e.
/usr/local/bin/hoppet-config). By default fastjet is not linked and only a
skeleton analysis (analysis/simple_analysis.f) is compiled. To link fastjet
run

⬇

cmake -DNEEDS_FASTJET=ON [-DFASTJET_CONFIG=/path/to/fastjet-config]

where the path to fastjet-config only needs to be specified if it is not in
the user’s $PATH. To compile a different analysis the user should first put it
in the analysis directory (here we assume it to be called my_analysis.f), and
then pass it to cmake through

⬇

cmake -DANALYSIS=my_analysis.f

The program has been found to compile on a Linux machine using the gfortran
v11.4.0 compiler and also on various MacOS systems.

####  3.1.1 The analysis framework

The code uses the POWHEG-BOX analysis framework, with some minor
modifications. A few example analyses are included in the analysis directory.
Any new analysis should be put here and the name of the analysis should be
passed to cmake as described above. There are two mandatory routines in the
analysis file, define_histograms and user_analysis. In the first routine one
should define histograms like this (there are also routines available to book
histograms with varying bins sizes)

⬇

call bookupeqbins(’string_name’, binsize, min, max)

The user_analysis routine takes as input

⬇

integer n

double precision dsig(maxscales), x, y, Q2

where n is the number of initial plus final state particles, dsig is the
weights computed by disorder and maxscales is the maximum number of scales
which is supported (currently 7). x, y, and Q2 are the DIS variables.

Through the module mod_analysis the analysis has access to two array of
momenta pbreit(0:3,1:n) and plab(0:3,1:n) in which the Breit and laboratory
frame momenta are stored respectively. The first entry is the incoming lepton,
the second the incoming parton, the third the outgoing lepton and the rest
outgoing partons. The output of the analysis will be saved to the disk as
outlined below. In the user_analysis routine the user should perform their
analysis and fill histograms like this

⬇

call filld(’string_name’, obs_value, dsig)

where obs_value is the value of the observable to be binned, and dsig is the
associated array of weights.

In the aux folder one may also find a script called mergedata that can perform
various manipulations of the datafiles. In particular

⬇

./mergedata 1 {list of statistically equivalent files}

will take the average of all files and produce the file fort.12 with the
result. Running the script without any arguments will result in a list of
possible uses of the script. The mergedata script is taken from the POWHEG-BOX
as well.

###  3.2 Inclusive mode

The syntax for running the program is

⬇

./disorder -pdf LHAPDF_name [options]

Running the program without any other options than -pdf will compute the total
inclusive cross section above the minimum $Q$ value accessible in the PDF, and
using default parameters everywhere. The user can get a list of most
parameters which can be specified on the command line by running

⬇

./disorder -help

Here we describe the most common flags, but for a complete list of all
parameters, and their use, the user should look through the file
src/mod_parameters.f90.

The program allows the user to specify limits on $x$, $y$, and $Q$ through
-xmin, -xmax, -ymin, -ymax, -Qmin, -Qmax, or to fix them through the options
-x, -y, -Q. For instance to compute the cross section at $Q=20\,\mathrm{GeV}$
and $x>0.01$ one would run

⬇

./disorder -Q 20 -xmin 0.01 -pdf LHAPDF_name

The program will perform a Monte Carlo integration in the ranges specified,
using the integrator VEGAS [64]. If the phase space is fully constrained by
fixing two of either $x$, $y$, and $Q$ the program simply evaluates one point
and returns the answer.

One can further specify the energy of the incoming lepton through -Elep, the
incoming hadron through -Ehad. By default the lepton is taken to be an
electron but specifying -positron on the command line will change that. The
code computes the photon-mediated NC cross section only by default. To include
the $Z$ one can specify -includeZ and to include CC processes one can specify
-CC. The inclusion of NC processes can also be controlled through the -NC
flag. In fact, all logical flags can be prefixed by “no” to turn them off.
Hence the below command line would run the program with CC processes only
using a positron

⬇

./disorder -Q 20 -xmin 0.01 -noNC -CC -positron -pdf LHAPDF_name

The order of the calculation is by default NNLO but can be specified with one
of the flags -lo/-nlo/-nnlo/-n3lo. If the user wants to compute PDF
uncertainties, the flag -pdfuncert should be given. This flag will make
disorder loop over all the members in the PDF, and combine their errors
according to the routine getpdfuncertainty [65], native to LHAPDF. This also
means that the program is slowed down proportionally to the number of PDF
members (although the VEGAS grid only gets computed once and then stored so
that the runs are fully correlated). If the PDF also includes
$\alpha_{\mathrm{s}}$ variations these are included in the PDF uncertainty by
default. If the user wants the PDF and $\alpha_{\mathrm{s}}$ uncertainties
independently, then the flag -alphasuncert should be specified. Some care
should be taken here, as in practice the code assumes that the
$\alpha_{\mathrm{s}}$ variations are contained in the last two PDF members,
and simply separates them from the rest.

Renormalisation, $\mu_{R}$, and factorisation, $\mu_{F}$, scale uncertainties
can be included by specifying the flag -scaleuncert. disorder uses the vector
boson virtuality, $Q^{2}$, as its default central scale, but can use any scale
as long as it is defined in terms of $x$, $y$, and $Q$. A few central scales,
$\mu$, are currently implemented and can be accessed through the -scale-choice
flag. The options are 1 (default): $Q^{2}$, 2: $Q^{2}(1-y)$, 3:
$\frac{Q^{2}(1-x)}{x}$. The program will then compute a standard 7-point scale
variation varying this scale by a factor of two up and down keeping
$1/2\leq\mu_{R}/\mu_{F}\leq 2$. On the screen the envelope of all 7 runs will
be printed. The user can in principle also carry out arbitrary variations in
individual runs by specifying -xmur and -xmuf on the command line. Here xmur
is the ratio of $\mu_{R}/\mu$ and similarly for xmuf.999In practice it is
faster to use the -scaleuncert flag as the program will only recompute what is
needed for the variations rather than do a full event. However, at N3LO, in
general, this is not true, due to the fact that the number of tables needed in
hoppet to carry out on-the-fly scale variations increases dramatically at this
order, compared to using a fixed ratio of $Q$. Given how fast the code is, it
is still the author’s opinion that it is more convenient to use the on-the-fly
variations.

Finally the random seed can be set with the iseed flag. When the program
terminates it will print results to screen but also save a number of files
depending on the exact input. There is always a file called
xsct_nnlo_seed0001.dat, where the nnlo and seed0001 parts will vary depending
on the order and seed, which contains a summary of the run, including the
total and reduced cross sections.

The output of the analysis is printed to a number of .dat files, the number
depending on the input, prefixed by disorder_nnlo_seed0001_pdfmem000 where
again the exact prefix will depend on the input. The output name also contains
information on the seed, the PDF member and the renormalisation and
factorisation scales if -scaleuncert is on. The user can specify a prefix to
be added to all the files through the -prefix flag.

Since it is often not necessary to run an analysis in the inclusive mode, the
user can simply turn this off with the -no-analysis flag.

Finally the user can control the number of VEGAS integration points through
the two flags -ncall1 and -ncall2. The first flags controls the number of
points to use to set up the grids. Since the phase space is not complicated it
is rarely necessary to increase the default of 10000. The second flag controls
the number of points that are used for the actual integration. If one is not
running an analysis, the default number of 100000 should give results that
have better than permille level accuracy. If the user wishes for instance to
bin the cross section in fine bins of $x$ and $Q^{2}$, this number will most
likely have to be increased.

###  3.3 P2B mode

To turn on P2B it is enough to specify -p2b on the command line. Many of the
flags described above can also be specified in this mode, with a few
limitations

  * •

-CC and -includeZ are not supported

  * •

-n3lo is not supported

  * •

-pdfuncert and -alphasuncert are not supported

Importantly, it _is_ possible to run with the -scaleuncert flag which leads to
a significant reduction in run-time compared to doing seven separate runs.
When running in P2B mode it is furthermore of use to be able to control the
number of calls to disent. This is done through the flag -ncall2, introduced
already above.

We provide a small script to run on multiple cores on a single machine in
aux/runpar.sh.

###  3.4 Validating the code

The code comes with a script to validate that the results come out as
expected. It can be run by entering validation and executing

⬇

./validate_or_generate.sh validate

The script builds the code and executes a number of tests designed to check
the most important features of disorder, including both the NC and CC
channels, the fastjet and LHAPDF interfaces, the automated PDF and scale
uncertainty features, and P2B. The script assumes default installation paths
for hoppet, fastjet, and LHAPDF. A user can however manually change this by
inspecting the script and adding the appropriate cmake flags at the beginning.
The validation script runs in about 500 CPU seconds, or no more than a few
minutes on a modern laptop.

##  4 Benchmarks and results

x

$Q=10\,\mathrm{GeV},\quad x=0.01$

|  $\frac{\sigma_{\mathrm{NC}}}{\mathrm{d}x\mathrm{d}Q^{2}}$ [pb/GeV2] |  $\delta$PDF |  $\frac{\sigma_{\mathrm{CC}}}{\mathrm{d}x\mathrm{d}Q^{2}}$ [pb/GeV2] |  $\delta$PDF | $\tilde{\sigma}_{\mathrm{NC}}$ | $\tilde{\sigma}_{\mathrm{CC}}$
---|---|---|---|---|---|---
aN3LO | $1932^{+.705\%}_{-.486\%}$ | $0.972\%$ | $1.100^{+.591\%}_{-.423\%}$ | $1.103\%$ | $0.8174$ | $2.884$
N3LO | $1886^{+.546\%}_{-.160\%}$ | $0.845\%$ | $1.080^{+.455\%}_{-.132\%}$ | $0.914\%$ | $0.7980$ | $2.829$
NNLO | $1895^{+1.50\%}_{-1.18\%}$ | $0.840\%$ | $1.084^{+1.25\%}_{-.994\%}$ | $0.908\%$ | $0.8018$ | $2.840$
NLO | $1952^{+3.66\%}_{-4.43\%}$ | $0.810\%$ | $1.111^{+3.02\%}_{-3.63\%}$ | $0.880\%$ | $0.8260$ | $2.913$
LO | $2058^{+13.8\%}_{-17.0\%}$ | $0.843\%$ | $1.163^{+11.4\%}_{-14.1\%}$ | $0.885\%$ | $0.9708$ | $3.050$

Table 1: The inclusive and reduced cross sections at various orders in both NC
and CC DIS. The setup is given in eq. (12). The aN3LO row is obtained using
the approximate N3LO PDF set MSHT20an3lo_as118. Note that there is no Monte
Carlo error on these numbers.

In this section we show a few results obtained by running the code in both the
inclusive and the P2B modes. The purpose is not to provide an exhaustive
phenomenological analysis, but rather show the capabilities of disorder and to
provide a few select results that can be used to either validate the code when
running it, or to validate other programs in the future. All analyses used in
this section can be found in the analysis directory, and the raw results can
be found in the paper_runs folder.

We use fastjet v3.4.1 for jet-clustering, and unless otherwise stated we
collide electrons and protons using the MSHT20nnlo_as118 [66] PDF set using
LHAPDF v6.5.4 and the following input parameters in all calculations

$\displaystyle M_{W}=80.398\,\mathrm{GeV},\qquad M_{Z}=91.1876\,\mathrm{GeV},$
$\displaystyle\qquad E_{h}=920\,\mathrm{GeV},\qquad E_{l}=27.6\,\mathrm{GeV},$
(12) $\displaystyle\alpha=1/137,\qquad\alpha_{\mathrm{s}}(M_{Z})$
$\displaystyle=0.118,\qquad n_{f}=5\,.$

###  4.1 Inclusive results

The perhaps most fundamental question one can ask in DIS is “what is the total
cross section for a given value of $x$ and $Q$?”. The answer to that question
is given in table 1 for both the NC and CC channels at $Q=10\,\mathrm{GeV}$
and $x=0.01$, at all available perturbative orders. The row aN3LO is obtained
with the MSHT20an3lo_as118 [67] PDF set which includes approximate N3LO
theoretical input. It is interesting to note that the inclusion of the
approximate terms has an impact that is parametrically of the same order as
the NNLO corrections. The PDF uncertainties stay more or less constant across
all orders, which is not unexpected given that we use the same PDF at each
order. To give a sense of the speed of the code, the combined results
presented in the table, which corresponds to 325 different PDF members and
hence structure functions, took a total of 3.5 minutes to obtain on a laptop
equipped with an Intel i9-10885H CPU.

Figure 1: The NC cross section for fixed $Q=10\,\mathrm{GeV}$ differential in
$\log x$ (left) and for fixed $x=0.01$ and differential in $\log Q$ (right).
We show NNLO (purple), N3LO (green), and aN3LO (red). The uncertainty band is
obtained as the linear combination of scale and PDF uncertainties.

In figure 1 we show the NC cross section for fixed $Q=10\,\mathrm{GeV}$ (left)
or $x=0.01$ (right), now only at NNLO, N3LO and aN3LO. The uncertainty bands
are obtained as a linear combination of the scale uncertainty and the PDF
uncertainty, and is typically dominated by the latter, in particular at N3LO.
We observe that although the central prediction at aN3LO is in tension with
the N3LO curve, it is contained within the NNLO scale uncertainty band, except
at intermediate $x$ values. It will be interesting to see, if the inclusion of
the exact N3LO splitting functions, when they become available, will
ameliorate this tension.

###  4.2 Exclusive laboratory frame results

To demonstrate the results of the code in P2B mode, we perform a jet analysis
in laboratory frame, based on Ref. [60], which probes the kinematic range
accessible by the upcoming EIC. For this analysis we set
$E_{l}=18\,\mathrm{GeV}$, $E_{h}=275\,\mathrm{GeV}$, and restrict the DIS
kinematics by

$\displaystyle 25\,\mathrm{GeV}^{2}<Q^{2}<1000\,\mathrm{GeV}^{2},\quad
0.04<y<0.95\,.$ (13)

We reconstruct jets in the laboratory frame using the anti-$k_{T}$ algorithm
[68] with $R=0.8$. Jets are those that satisfy the following transverse
momentum and pseudo-rapidity requirements

$\displaystyle p_{t,j}>5\,\mathrm{GeV},\quad|\eta_{j}|<3\,.$ (14)

Figure 2: The NC cross section for fixed $Q=10\,\mathrm{GeV}$ differential in
$\log x$ (left) and for fixed $x=0.01$ and differential in $\log Q$ (right).
We show NNLO (purple), N3LO (green), and aN3LO (red). The uncertainty band is
obtained as the linear combination of scale and PDF uncertainties.

In Fig. 2 we show the inclusive jet transverse momentum and rapidity at LO,
NLO, and NNLO, in photon-induced NC DIS. For the transverse momentum the
higher order corrections are moderate except for small values and close to the
LO threshold of $p_{t,j}=Q$. The rapidity distribution on the other hand
receives larger corrections. This is in part due to hard jets with $p_{t,j}>Q$
that can only arise at NLO and beyond. Here the scale variation bands do not
cover the actual corrections, and it is obvious that one needs to go beyond
NLO for reliable predictions. Of course one needs the full N3LO to assess the
reliability of the NNLO result. However, from similar results in Ref. [36] it
can be expected that the actual N3LO results will be mostly contained within
the scale variation band at NNLO.

##  5 Conclusion

In this paper we have presented version 1.0.0 of a Fortran code, disorder,
capable of computing massless DIS cross sections fully differentially at NNLO
and inclusively at N3LO. The code combines the structure functions from hoppet
and the exclusive DIS + 1 jet calculation of disent, using the P2B method. In
the inclusive mode the program includes all NC and CC contributions and can
carry out both scale and PDF variations on the fly. Due to the underlying
tabulation of the structure functions in hoppet, the code is extremely fast,
and a user can obtain cross section results at N3LO in a matter of seconds.

In the fully differential P2B mode, only photon-mediated DIS can be computed,
and no PDF variations are currently implemented. It should in principle be
possible to extend disent to full NC and CC, and this is something that is
planned for a future release of the code. Similarly the PDF variations could
also be implemented, however this would require some serious restructuring of
the code for it to be efficient. Extension of the code to handle incoming
neutrino beams, of relevance to the FASER [69] physics program, is also
planned.

The upgrade of the code to handle N3LO fully differentially is however not
currently planned. This is mainly due to the fact that this would require a
complete replacement of disent with a proper exclusive NNLO DIS + 1 jet code,
entailing a significant amount of work. Such a code would also be very slow
compared to what is in place now, and it would defeat part of the attraction
of disorder which is its speed. At LO and NLO a user can obtain differential
results on a laptop running for a just a few minutes, and even reasonable NNLO
results can be obtained on a multi-core machine running for an hour. N3LO
results would inevitably require large resources on a High Performance
Cluster.

Given the code’s dependency on hoppet, any developments therein will almost
automatically propagate to disorder. In particular it is expected that the
exact $\mathcal{O}(\alpha_{\mathrm{s}}^{4})$ splitting functions, as opposed
to the approximations which are currently employed, will become available in
hoppet once they have been fully determined. This is currently the only
missing piece needed to formally claim full N3LO accuracy for the structure
functions. Similarly, were the massive DIS coefficient functions [70, 71, 72,
73, 74, 75, 76, 77] to be implemented in hoppet the structure functions in
disorder could immediately be modified to accommodate this.

The code is intended to be user-friendly and to this effect comes with an
interface to fastjet and LHAPDF and is run through a command line interface.

## Acknowledgements

The author is grateful for many discussions with Andrea Banfi and Gavin Salam
regarding disent, and for encouragement to publish disorder. The author also
thanks Valerio Bertone for discussions about the DIS structure functions and
Alex Huss for valuable cross-checks against results obtained with the NNLOjet
code. The idea for the code originated while the author was working on Ref.
[78].

## References

  * [1] J. C. Collins and D. E. Soper,  _The Theorems of Perturbative QCD_ ,  Ann. Rev. Nucl. Part. Sci. 37, 383 (1987),  10.1146/annurev.ns.37.120187.002123.
  * [2] J. C. Collins, D. E. Soper and G. F. Sterman,  _Factorization of Hard Processes in QCD_ ,  Adv. Ser. Direct. High Energy Phys. 5, 1 (1989),  10.1142/9789814503266_0001,  hep-ph/0409313.
  * [3] J. Sanchez Guillen, J. Miramontes, M. Miramontes, G. Parente and O. A. Sampayo,  _Next-to-leading order analysis of the deep inelastic R = sigma-L / sigma-total_ ,  Nucl. Phys. B 353, 337 (1991),  10.1016/0550-3213(91)90340-4.
  * [4] W. L. van Neerven and E. B. Zijlstra,  _Order alpha-s**2 contributions to the deep inelastic Wilson coefficient_ ,  Phys. Lett. B 272, 127 (1991),  10.1016/0370-2693(91)91024-P.
  * [5] E. B. Zijlstra and W. L. van Neerven,  _Order alpha-s**2 QCD corrections to the deep inelastic proton structure functions F2 and F(L)_ ,  Nucl. Phys. B 383, 525 (1992),  10.1016/0550-3213(92)90087-R.
  * [6] E. B. Zijlstra and W. L. van Neerven,  _Order alpha-s**2 correction to the structure function F3 (x, Q**2) in deep inelastic neutrino - hadron scattering_ ,  Phys. Lett. B 297, 377 (1992),  10.1016/0370-2693(92)91277-G.
  * [7] W. L. van Neerven and A. Vogt,  _NNLO evolution of deep inelastic structure functions: The Nonsinglet case_ ,  Nucl. Phys. B 568, 263 (2000),  10.1016/S0550-3213(99)00668-9,  hep-ph/9907472.
  * [8] W. L. van Neerven and A. Vogt,  _NNLO evolution of deep inelastic structure functions: The Singlet case_ ,  Nucl. Phys. B 588, 345 (2000),  10.1016/S0550-3213(00)00480-6,  hep-ph/0006154.
  * [9] S. Moch and J. A. M. Vermaseren,  _Deep inelastic structure functions at two loops_ ,  Nucl. Phys. B 573, 853 (2000),  10.1016/S0550-3213(00)00045-6,  hep-ph/9912355.
  * [10] S. Moch, J. A. M. Vermaseren and A. Vogt,  _The Longitudinal structure function at the third order_ ,  Phys. Lett. B 606, 123 (2005),  10.1016/j.physletb.2004.11.063,  hep-ph/0411112.
  * [11] J. A. M. Vermaseren, A. Vogt and S. Moch,  _The Third-order QCD corrections to deep-inelastic scattering by photon exchange_ ,  Nucl. Phys. B 724, 3 (2005),  10.1016/j.nuclphysb.2005.06.020,  hep-ph/0504242.
  * [12] A. Vogt, S. Moch and J. Vermaseren,  _Third-order QCD results on form factors and coefficient functions_ ,  Nucl. Phys. B Proc. Suppl. 160, 44 (2006),  10.1016/j.nuclphysbps.2006.09.101,  hep-ph/0608307.
  * [13] S. Moch, M. Rogal and A. Vogt,  _Differences between charged-current coefficient functions_ ,  Nucl. Phys. B 790, 317 (2008),  10.1016/j.nuclphysb.2007.09.022,  0708.3731.
  * [14] J. Davies, A. Vogt, S. Moch and J. A. M. Vermaseren,  _Non-singlet coefficient functions for charged-current deep-inelastic scattering to the third order in QCD_ ,  PoS DIS2016, 059 (2016),  10.22323/1.265.0059,  1606.08907.
  * [15] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald,  _The massless three-loop Wilson coefficients for the deep-inelastic structure functions F 2, FL, xF3 and g1_,  JHEP 11, 156 (2022),  10.1007/JHEP11(2022)156,  2208.14325.
  * [16] S.-O. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,  _DIS coefficient functions at four loops in QCD and beyond_ ,  PoS LL2022, 047 (2022),  10.22323/1.416.0047,  2208.11067.
  * [17] S. Moch, J. A. M. Vermaseren and A. Vogt,  _The Three loop splitting functions in QCD: The Nonsinglet case_ ,  Nucl. Phys. B 688, 101 (2004),  10.1016/j.nuclphysb.2004.03.030,  hep-ph/0403192.
  * [18] A. Vogt, S. Moch and J. A. M. Vermaseren,  _The Three-loop splitting functions in QCD: The Singlet case_ ,  Nucl. Phys. B 691, 129 (2004),  10.1016/j.nuclphysb.2004.04.024,  hep-ph/0404111.
  * [19] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald,  _The three-loop unpolarized and polarized non-singlet anomalous dimensions from off shell operator matrix elements_ ,  Nucl. Phys. B 971, 115542 (2021),  10.1016/j.nuclphysb.2021.115542,  2107.06267.
  * [20] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin,  _The Four loop beta function in quantum chromodynamics_ ,  Phys. Lett. B 400, 379 (1997),  10.1016/S0370-2693(97)00370-5,  hep-ph/9701390.
  * [21] M. Czakon,  _The Four-loop QCD beta-function and anomalous dimensions_ ,  Nucl. Phys. B 710, 485 (2005),  10.1016/j.nuclphysb.2005.01.012,  hep-ph/0411261.
  * [22] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,  _Low moments of the four-loop splitting functions in QCD_ ,  Phys. Lett. B 825, 136853 (2022),  10.1016/j.physletb.2021.136853,  2111.15561.
  * [23] G. Falcioni, F. Herzog, S. Moch and A. Vogt,  _Four-loop splitting functions in QCD – The quark-quark case_ ,  Phys. Lett. B 842, 137944 (2023),  10.1016/j.physletb.2023.137944,  2302.07593.
  * [24] G. Falcioni, F. Herzog, S. Moch and A. Vogt,  _Four-loop splitting functions in QCD – The gluon-to-quark case_ (2023),  2307.04158.
  * [25] T. Gehrmann, A. von Manteuffel, V. Sotnikov and T.-Z. Yang,  _Complete $N_{f}^{2}$ contributions to four-loop pure-singlet splitting functions_ (2023),  2308.07958.
  * [26] G. Falcioni, F. Herzog, S. Moch, J. Vermaseren and A. Vogt,  _The double fermionic contribution to the four-loop quark-to-gluon splitting function_ (2023),  2310.01245.
  * [27] S. Moch, B. Ruijl, T. Ueda, J. Vermaseren and A. Vogt,  _Additional moments and x-space approximations of four-loop splitting functions in QCD_ (2023),  2310.05744.
  * [28] T. Gehrmann, A. von Manteuffel, V. Sotnikov and T.-Z. Yang,  _The $N_{f}\,C_{F}^{3}$ contribution to the non-singlet splitting function at four-loop order_ (2023),  2310.12240.
  * [29] I. Bierenbaum, J. Blumlein and S. Klein,  _Mellin Moments of the O(alpha**3(s)) Heavy Flavor Contributions to unpolarized Deep-Inelastic Scattering at Q**2 $>$$>$ m**2 and Anomalous Dimensions_,  Nucl. Phys. B 820, 417 (2009),  10.1016/j.nuclphysb.2009.06.005,  0904.3563.
  * [30] H. Kawamura, N. A. Lo Presti, S. Moch and A. Vogt,  _On the next-to-next-to-leading order QCD corrections to heavy-quark production in deep-inelastic scattering_ ,  Nucl. Phys. B 864, 399 (2012),  10.1016/j.nuclphysb.2012.07.001,  1205.5727.
  * [31] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider and F. Wißbrock,  _The transition matrix element $a_{gq}(n)$ of the variable flavor number scheme at $\mathcal{O}(\alpha_{s}^{3})$_,  Nuclear Physics B 882, 263 (2014),  https://doi.org/10.1016/j.nuclphysb.2014.02.007,  1402.0359.
  * [32] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider and F. Wißbrock,  _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_,  Nucl. Phys. B 886, 733 (2014),  10.1016/j.nuclphysb.2014.07.010,  1406.4654.
  * [33] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider,  _The 3-loop pure singlet heavy flavor contributions to the structure function $F_{2}(x,Q^{2})$ and the anomalous dimension_,  Nucl. Phys. B 890, 48 (2014),  10.1016/j.nuclphysb.2014.10.008,  1409.1135.
  * [34] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Goedicke, A. von Manteuffel, C. Schneider and K. Schönwald,  _The unpolarized and polarized single-mass three-loop heavy flavor operator matrix elements A gg,Q and $\Delta$Agg,Q_,  JHEP 12, 134 (2022),  10.1007/JHEP12(2022)134,  2211.05462.
  * [35] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider and K. Schönwald,  _The first-order factorizable contributions to the three-loop massive operator matrix elements $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$_ (2023),  2311.00644.
  * [36] J. Currie, T. Gehrmann, E. W. N. Glover, A. Huss, J. Niehues and A. Vogt,  _N 3LO corrections to jet production in deep inelastic scattering using the Projection-to-Born method_,  JHEP 05, 209 (2018),  10.1007/JHEP05(2018)209,  1803.09973.
  * [37] T. Gehrmann, A. Huss, J. Niehues, A. Vogt and D. M. Walker,  _Jet production in charged-current deep-inelastic scattering to third order in QCD_ ,  Phys. Lett. B 792, 182 (2019),  10.1016/j.physletb.2019.03.003,  1812.06104.
  * [38] M. Cacciari, F. A. Dreyer, A. Karlberg, G. P. Salam and G. Zanderighi,  _Fully Differential Vector-Boson-Fusion Higgs Production at Next-to-Next-to-Leading Order_ ,  Phys. Rev. Lett. 115(8), 082002 (2015),  10.1103/PhysRevLett.115.082002,  [Erratum: Phys.Rev.Lett. 120, 139901 (2018)],  1506.02660.
  * [39] E. Mirkes and D. Zeppenfeld,  _Dijet production at HERA in next-to-leading order_ ,  Phys. Lett. B 380, 205 (1996),  10.1016/0370-2693(96)00426-1,  hep-ph/9511448.
  * [40] S. Catani and M. H. Seymour,  _A General algorithm for calculating jet cross-sections in NLO QCD_ ,  Nucl. Phys. B 485, 291 (1997),  10.1016/S0550-3213(96)00589-5,  [Erratum: Nucl.Phys.B 510, 503–504 (1998)],  hep-ph/9605323.
  * [41] D. Graudenz,  _Disaster++: Version 1.0_ (1997),  hep-ph/9710244.
  * [42] Z. Nagy and Z. Trocsanyi,  _Multijet cross-sections in deep inelastic scattering at next-to-leading order_ ,  Phys. Rev. Lett. 87, 082001 (2001),  10.1103/PhysRevLett.87.082001,  hep-ph/0104315.
  * [43] V. Antonelli, M. Dasgupta and G. P. Salam,  _Resummation of thrust distributions in DIS_ ,  JHEP 02, 001 (2000),  10.1088/1126-6708/2000/02/001,  hep-ph/9912488.
  * [44] M. Dasgupta and G. P. Salam,  _Resummed event shape variables in DIS_ ,  JHEP 08, 032 (2002),  10.1088/1126-6708/2002/08/032,  hep-ph/0208073.
  * [45] I. Borsa, D. de Florian and I. Pedron,  _Jet Production in Polarized Deep Inelastic Scattering at Next-to-Next-to-Leading Order_ ,  Phys. Rev. Lett. 125(8), 082001 (2020),  10.1103/PhysRevLett.125.082001,  2005.10705.
  * [46] I. Borsa, D. de Florian and I. Pedron,  _Inclusive-jet and dijet production in polarized deep inelastic scattering_ ,  Phys. Rev. D 103(1), 014008 (2021),  10.1103/PhysRevD.103.014008,  2010.07354.
  * [47] G. P. Salam and J. Rojo,  _A Higher Order Perturbative Parton Evolution Toolkit (HOPPET)_ ,  Comput. Phys. Commun. 180, 120 (2009),  10.1016/j.cpc.2008.08.010,  0804.3755.
  * [48] V. Bertone and A. Karlberg,  _Benchmark of the deep-inelastic-scattering structure functions at $\mathcal{O}(\alpha_{\mathrm{s}}^{3})$_,  CERN-TH-2024-018.
  * [49] V. Bertone, S. Carrazza and J. Rojo,  _APFEL: A PDF Evolution Library with QED corrections_ ,  Comput. Phys. Commun. 185, 1647 (2014),  10.1016/j.cpc.2014.03.007,  1310.1394.
  * [50] V. Bertone,  _APFEL++: A new PDF evolution library in C++_ ,  PoS DIS2017, 201 (2018),  10.22323/1.297.0201,  1708.00911.
  * [51] F. A. Dreyer and A. Karlberg,  _Vector-Boson Fusion Higgs Production at Three Loops in QCD_ ,  Phys. Rev. Lett. 117(7), 072001 (2016),  10.1103/PhysRevLett.117.072001,  1606.00840.
  * [52] F. A. Dreyer and A. Karlberg,  _Vector-Boson Fusion Higgs Pair Production at N 3LO_,  Phys. Rev. D 98(11), 114016 (2018),  10.1103/PhysRevD.98.114016,  1811.07906.
  * [53] F. A. Dreyer and A. Karlberg,  _Fully differential Vector-Boson Fusion Higgs Pair Production at Next-to-Next-to-Leading Order_ ,  Phys. Rev. D 99(7), 074028 (2019),  10.1103/PhysRevD.99.074028,  1811.07918.
  * [54] A. Karlberg, P. Nason, G. Salam, G. Zanderighi and F. Dreyer,  _Hoppet v1.3.0 release note_ ,  CERN-TH-2023-237, MPP-2023-285, OUTP-23-15P.
  * [55] M. Cacciari, G. P. Salam and G. Soyez,  _FastJet User Manual_ ,  Eur. Phys. J. C 72, 1896 (2012),  10.1140/epjc/s10052-012-1896-2,  1111.6097.
  * [56] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt,  _LHAPDF6: parton density access in the LHC precision era_ ,  Eur. Phys. J. C 75, 132 (2015),  10.1140/epjc/s10052-015-3318-8,  1412.7420.
  * [57] 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, 043 (2010),  10.1007/JHEP06(2010)043,  1002.2581.
  * [58] I. Borsa, D. de Florian and I. Pedron,  _Dijet production in neutral current and charged current polarized deep inelastic scattering_ ,  Phys. Rev. D 105(7), 074025 (2022),  10.1103/PhysRevD.105.074025,  2112.08223.
  * [59] I. Borsa, D. de Florian and I. Pedron,  _The full set of polarized deep inelastic scattering structure functions at NNLO accuracy_ ,  Eur. Phys. J. C 82(12), 1167 (2022),  10.1140/epjc/s10052-022-11140-z,  2210.12014.
  * [60] I. Borsa, D. de Florian and I. Pedron,  _NNLO jet production in neutral and charged current polarized deep inelastic scattering_ ,  Phys. Rev. D 107(5), 054027 (2023),  10.1103/PhysRevD.107.054027,  2212.06625.
  * [61] R. Abdul Khalek _et al._ ,  _Snowmass 2021 White Paper: Electron Ion Collider for High Energy Physics_ (2022),  2203.13199.
  * [62] R. Devenish and A. Cooper-Sarkar,  _Deep inelastic scattering_ ,  10.1093/acprof:oso/9780198506713.001.0001 (2004).
  * [63] F. D. Aaron _et al._ ,  _Inclusive Deep Inelastic Scattering at High $Q^{2}$ with Longitudinally Polarised Lepton Beams at HERA_,  JHEP 09, 061 (2012),  10.1007/JHEP09(2012)061,  1206.7007.
  * [64] G. P. Lepage,  _A New Algorithm for Adaptive Multidimensional Integration_ ,  J. Comput. Phys. 27, 192 (1978),  10.1016/0021-9991(78)90004-9.
  * [65] G. Watt,  _Parton distribution function dependence of benchmark Standard Model total cross sections at the 7 TeV LHC_ ,  JHEP 09, 069 (2011),  10.1007/JHEP09(2011)069,  1106.5788.
  * [66] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne,  _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ ,  Eur. Phys. J. C 81(4), 341 (2021),  10.1140/epjc/s10052-021-09057-0,  2012.04684.
  * [67] J. McGowan, T. Cridge, L. A. Harland-Lang and R. S. Thorne,  _Approximate N 3LO parton distribution functions with theoretical uncertainties: MSHT20aN3LO PDFs_,  Eur. Phys. J. C 83(3), 185 (2023),  10.1140/epjc/s10052-023-11236-0,  [Erratum: Eur.Phys.J.C 83, 302 (2023)],  2207.04739.
  * [68] M. Cacciari, G. P. Salam and G. Soyez,  _The anti- $k_{t}$ jet clustering algorithm_,  JHEP 04, 063 (2008),  10.1088/1126-6708/2008/04/063,  0802.1189.
  * [69] J. L. Feng, I. Galon, F. Kling and S. Trojanowski,  _ForwArd Search ExpeRiment at the LHC_ ,  Phys. Rev. D 97(3), 035001 (2018),  10.1103/PhysRevD.97.035001,  1708.09389.
  * [70] T. Gottschalk,  _Chromodynamic Corrections to Neutrino Production of Heavy Quarks_ ,  Phys. Rev. D 23, 56 (1981),  10.1103/PhysRevD.23.56.
  * [71] E. Laenen, S. Riemersma, J. Smith and W. L. van Neerven,  _Complete O (alpha-s) corrections to heavy flavor structure functions in electroproduction_ ,  Nucl. Phys. B 392, 162 (1993),  10.1016/0550-3213(93)90201-Y.
  * [72] E. Laenen, S. Riemersma, J. Smith and W. L. van Neerven,  _O(alpha-s) corrections to heavy flavor inclusive distributions in electroproduction_ ,  Nucl. Phys. B 392, 229 (1993),  10.1016/0550-3213(93)90202-Z.
  * [73] M. Gluck, S. Kretzer and E. Reya,  _Detailed next-to-leading order analysis of deep inelastic neutrino induced charm production off strange sea partons_ ,  Phys. Lett. B 398, 381 (1997),  10.1016/S0370-2693(97)00232-3,  [Erratum: Phys.Lett.B 405, 392 (1997)],  hep-ph/9701364.
  * [74] J. Blumlein, A. Hasselhuhn, P. Kovacikova and S. Moch,  _$O(\alpha_{s})$ Heavy Flavor Corrections to Charged Current Deep-Inelastic Scattering in Mellin Space_,  Phys. Lett. B 700, 294 (2011),  10.1016/j.physletb.2011.05.007,  1104.3449.
  * [75] A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel and C. Schneider,  _O( $\alpha_{s}^{3}$) heavy flavor contributions to the charged current structure function $xF_{3}(x,Q^{2})$ at large momentum transfer_,  Phys. Rev. D 92(11), 114005 (2015),  10.1103/PhysRevD.92.114005,  1508.01449.
  * [76] E. L. Berger, J. Gao, C. S. Li, Z. L. Liu and H. X. Zhu,  _Charm-Quark Production in Deep-Inelastic Neutrino Scattering at Next-to-Next-to-Leading Order in QCD_ ,  Phys. Rev. Lett. 116(21), 212002 (2016),  10.1103/PhysRevLett.116.212002,  1601.05430.
  * [77] J. Gao,  _Massive charged-current coefficient functions in deep-inelastic scattering at NNLO and impact on strange-quark distributions_ ,  JHEP 02, 026 (2018),  10.1007/JHEP02(2018)026,  1710.04258.
  * [78] A. Banfi, S. F. Ravasio, B. Jäger, A. Karlberg, F. Reichenbach and G. Zanderighi,  _A POWHEG generator for deep inelastic scattering_ (2023),  2309.02127.