# The role of multi-parton interactions in doubly-heavy hadron production

U. Egede1, T. Hadavizadeh1, M. Singla1, P. Skands1,
M. Vesterinen2
1 _School of Physics and Astronomy, Monash University, Melbourne, Australia_
2 _Department of Physics, University of Warwick, Coventry, United Kingdom_

(August 2022)

###### Abstract

Beauty and charm quarks are ideal probes of pertubative Quantum Chromodymanics
in proton-proton collisions, owing to their large masses. In this paper the
role of multi-parton interactions in the production of doubly-heavy hadrons is
studied using simulation samples generated with Pythia, a Monte Carlo event
generator. Comparisons are made to the stand-alone generators Bcvegpy and
GenXicc. New methods of speeding up Pythia simulations for events containing
heavy quarks are described, enabling the production of large samples with
multiple heavy-quark pairs. We show that significantly higher production rates
of doubly-heavy hadrons are predicted in models that allow heavy quarks
originating from different parton-parton interactions (within the same hadron-
hadron collision) to combine to form such hadrons. Quantitative predictions
are sensitive to the modelling of colour reconnections. We suggest a set of
experimental measurements capable of differentiating these additional
contributions.

© The authors, CC-BY-4.0 licence.

##  1 Introduction

The masses of the charm and beauty quarks are so high that their production at
hadron colliders is dominated by perturbative QCD processes making them ideal
probes of the partonic interactions that led to their formation. Hadrons
containing two beauty and/or charm quarks are referred to as doubly-heavy
hadrons and, here, include both mesons and baryons with a net heavy flavour
such as ${B}_{c}^{+}$ mesons and ${\mathchar 28932\relax}^{++}_{{c}{c}}$
baryons, and quarkonia with zero net flavour.

Large samples of doubly-heavy hadrons have been collected at the Large Hadron
Collider (LHC) and the study of their properties can provide a unique insight
into the role of multi-parton interactions (MPIs) in hadron formation. Because
$c$ and $b$ quarks are too heavy to be produced non-perturbatively, doubly-
heavy hadrons can only be formed by the effective coalescence of two
perturbatively produced heavy quarks. This is fundamentally different from
singly-heavy hadrons, in which the heavy quark is confined together with light
quark(s) that can be produced non-perturbatively. The fact that heavy hadrons
can be produced via MPI has been confirmed not only in inclusive measurements
of heavy-hadron-pair cross sections [1, 2, 3, 4] but also in differential
measurements of ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ rates vs
charged-track multiplicity [5]. However, these measurements do not by
themselves show conclusively whether two heavy quarks from _different_ parton-
parton interactions can join to form hadrons.111Throughout, all parton-parton
interactions are within the context of a single proton-proton collision.
Recent measurements of the newly-discovered doubly-charmed tetraquark
$T_{cc}^{++}$ indicate the production has similarities with processes
involving different parton-parton interactions [6]. In this paper, we point
out that this question can be addressed by considering non-onium doubly-heavy
hadrons, such as $B_{c}^{+}$ and $\Xi_{cc}^{++}$. Moreover, such measurements
will place constraints on models of colour reconnections (CR), which are
relevant to a broad range of hadron-collider physics studies [7, 8, 9, 10].

Samples of doubly-heavy hadrons can be simulated using Monte Carlo event
generators such as Pythia [11, 12]. However, as these hadrons require the
chance coalescence of two heavy particles into a bound state, their generation
can be prohibitively slow. For this reason some rare doubly-heavy hadrons, for
example, ${B}_{c}^{+}$, $\Xi_{cc}$ and $\Xi_{bc}$, are at the moment usually
generated using dedicated generators, such as Bcvegpy [13] or GenXicc [14],
that perform fixed-order matrix-element calculations. These generators are
then interfaced with event generators to simulate the rest of the event
evolution and hadronisation. However, the fixed-order matrix-element
generators assume that both of the heavy quarks are produced in a _single_
parton-parton interaction, for example $gg\rightarrow B_{c}^{+}b\bar{c}$ or
$q\bar{q}\rightarrow B_{c}^{+}b\bar{c}$. This ignores the role that MPIs could
play in generating heavy quarks that could contribute to the formation of such
hadrons.

In this paper, we develop a method to enhance the rate of hadrons containing
one or more heavy quarks in simulation samples generated with Pythia. This
enables us to investigate the role of MPIs in the production of doubly-heavy
hadrons in the context of inclusive event samples, in which the formation of
such hadrons would normally be exceedingly rare. The generated samples are
studied to identify properties that can differentiate between the
contributions from single parton scattering (SPS) and double parton scatering
(DPS) mechanisms. Simulations are carried out with Pythia version 8.306 using
the default Monash tune [15] and at a proton-proton centre of mass energy of
$13\mathrm{\,Te\kern-1.00006ptV}$. The hard interaction is simulated using
Pythia 8.306, Bcvegpy 2.2 or GenXicc 2.1 while the fragmentation process is
simulated using the default Pythia Simple Shower framework [16, 17]. Finally,
measurements using LHC data that are able to shed light on the production
mechanism are proposed and their feasibility discussed.

Throughout this paper charge conjugation is not implied when describing cross
sections. The generator Bcvegpy only produces ${B}_{c}^{+}$ mesons, therefore
to avoid ambiguity all cross sections are explicitly only produced for the
explicitly specified charge. In general the cross section predictions from
Bcvegpy should be doubled when comparing to measurements of
$\sigma(\mbox{${p}{p}\\!\rightarrow B_{c}^{\pm}X$})$.

##  2 Sources of heavy quarks

(a) Pair creation.

(b) Flavour excitation.

(c) Gluon splitting in parton shower.

Figure 1: Examples of production mechanisms for heavy quarks in proton-proton
collisions. The incoming, outgoing and intermediate particles of the process
considered to be the hardest process are highlighted in red. In the case of
flavour excitation, the $\bar{b}$ quark shown at the bottom represents the
companion quark produced as a result of the initial-state evolution.

In proton-proton collisions, the QCD production mechanisms for the heavy $c$
and $b$ quarks can be split into three categories referred to as pair
creation, flavour excitation and gluon splitting in parton showers [18]. The
processes are classified according to the interaction with the largest
momentum transfer, referred to here as the hard interaction.

Pair creation involves a $gg\rightarrow Q\bar{Q}$222$Q$ here represents heavy
quarks and $q$ represents light quarks or $q\bar{q}\rightarrow Q\bar{Q}$ hard
interaction, as shown in Fig. 1(a), that, in the absence of significant
initial-state radiation, creates outgoing heavy quarks with equal and opposite
transverse momenta. The resulting heavy hadrons formed from the heavy quarks
similarly have a strong tendency to be back-to-back in the transverse plane,
as shown in Fig. 2 for $b\bar{b}$ production.

Figure 2: Differential cross-section in the transverse angle between the two
$b$-hadrons in Pythia simulations with exactly two $b$-hadrons.

Flavour excitation is the process involving one heavy quark: $Qg\rightarrow
Qg$ or $Qq\rightarrow Qq$, represented in Fig. 1(b). In this process a virtual
$Q\bar{Q}$ pair is produced as part of the initial-state evolution of one of
the incoming protons, and one of them, say the $Q$, subsequently interacts
with a (non-heavy) parton from the other proton. The $\bar{Q}$ (a.k.a. the
“companion” quark of the scattered heavy quark [19]) is ejected as part of the
initial-state evolution of the incoming remnant at a lower scale, with less
transverse momentum and significantly less correlation with the direction of
the $Q$, as shown in Fig. 2.

Heavy quarks can also be produced via gluon splittings during parton showers.
A typical example would be a hard $gg\rightarrow gg$ interaction followed by a
subsequent $g\rightarrow Q\bar{Q}$ splitting in the subsequent initial- or
final-state shower evolution, as shown in Fig. 1(c). Although this figure
shows one of the outgoing gluons from the hard interaction directly splitting
to heavy quarks, that is just for simplicity; in principle any gluon produced
within a shower above the heavy quark-mass threshold could result in heavy
quarks. As gluon-gluon interactions have a large cross-section at the LHC,
this constitutes a significant contribution to the heavy-quark production
mechanisms. For final-state gluon splittings, the resulting $Q\bar{Q}$ pair
will be boosted in the direction of the parent gluon. Events in which two
singly-heavy hadron are produced by this mechanism tend to have smaller angles
between the two heavy hadrons, as shown in Fig. 2.

##  3 Sources of doubly-heavy hadrons

To create doubly-heavy hadrons that are not quarkonium states, two $Q\bar{Q}$
pairs must be produced during the perturbative evolution of the collision. An
example of an SPS mechanism contributing to this process is shown in Fig.
3(a): hard $b\bar{b}$ pair creation followed by a $g\rightarrow c\bar{c}$
splitting during the shower evolution. Equivalent processes involving flavour
excitation or double gluon splitting within a single SPS are of course also
possible.

When allowing for MPI, the two $Q\bar{Q}$ pairs may also be produced in two
_different_ parton-parton interactions (still within the context of a single
hadron-hadron collision). This is what we label DPS. Two examples, double pair
creation and double flavour excitation, are shown in Figs. 3(b) and 3(c)
respectively, again with other combinations of pair creation, flavour
excitation, and/or gluon splittings obviously also possible. In these diagrams
the two parton interactions have been highlighted in different colours to
clarify the origin of the partons.

In events with more than two parton-parton interactions, SPS mechanisms could
contribute from any one of the _single_ parton-parton interactions, whilst DPS
mechanisms could contribute from the combination of any two.

(a) Example of SPS: Pair creation and gluon splitting.

(b) Example of DPS: Double pair creation.

(c) Example of DPS: Double flavour excitation.

Figure 3: Production mechanisms for events with both a $b\bar{b}$ and
$c\bar{c}$ pair. The incoming, outgoing and intermediate particles of each
parton-parton interaction are shown in red and (where relevant) blue. In the
case of double flavour excitation, $b$ and $\bar{c}$ quarks shown at the top
and bottom represents the companion quark produced as a result of the initial-
state evolution.

Once the appropriate quarks have been produced in the collision, only pairs
that are sufficiently close in phase space and which have a non-zero
probability to be in an overall colour-singlet state, have a chance to form an
on-shell doubly-heavy hadron.

Measurements of the cross sections of multiple heavy hadrons suggest that MPIs
play a significant role in the production of multiple heavy quark pairs at
hadron colliders [5, 1, 2, 3, 4]. However, the question of how partons
originating from _different_ parts of the protons become bound into hadrons is
still afflicted with significant uncertainties. In general-purpose event
generators like Pythia, this is controlled by a combination of perturbative
heavy-quark production mechanisms (hard scatterings, MPI, and parton showers)
and semi-empirical models of colour reconnections with [20, 10] and without
[7, 8] space-time dependence. The simple diagrams in Fig. 3 demonstrate how
$B_{c}^{+}$ mesons formed from the $\bar{b}c$ combinations could provide an
ideal probe into the hadronisation process. This is unique to doubly-heavy
hadrons, since light quarks are mainly created nonperturbatively and hence do
not have the same character of being associated with specific short-distance
processes in the colliding protons.

##  4 Efficient simulation of events with heavy hadrons in Pythia

Generating unbiased events with multiple pairs of heavy quarks and doubly-
heavy hadrons with Monte Carlo event generators can be very time consuming as
few events will fulfil the requirements to form the doubly-heavy hadrons. A
method of enhancing the efficiency to produce events containing heavy quarks
in Pythia is outlined here, and can be applied to both singly- and doubly-
heavy hadrons.

Pythia provides user-configurable classes called UserHooks aimed at allowing
the user to inspect and veto events at different stages during the event
evolution. These can be exploited to veto events that do not contain the
requisite heavy quarks early on in the generation, removing time spent
evolving and hadronising events that will never be accepted.

(a) Process Level.

(b) During evolution.

(c) Parton Level.

(d) Hadronisation.

Figure 4: Simple representations of different stages during the event
evolution in Pythia.

The UserHook stages that are utilised to improve the efficiency are:

  * •

Hard-process-level veto: This veto inspects the event after the most energetic
parton interaction has occurred, as represented in Fig 4(a);

  * •

Event-evolution veto: In Pythia the event is evolved from the hard-interaction
scale down to the hadronisation scale. During this process, the event can be
inspected when the evolution reaches an arbitrary (user-defined) value of the
evolution scale, illustrated in Fig 4(b);

  * •

Parton-level veto: Once the evolution has reached the hadronisation scale, the
event can be inspected before any hadrons are created, leaving the complete
partonic event as shown in Fig 4(c);

  * •

Post-hadronisation veto: It is possible to inspect the event after the
hadronic states have been formed as shown in Fig 4(d).

Obviously, vetoes that can be placed early on can lead to far bigger
efficiency gains than ones placed later, and a veto at the post-hadronisation
level, after most of the event activity has already been treated (apart from
decays and, optionally, hadronic rescatterings) will not make much of a
difference at all.

###  4.1 Simulating final states involving one $Q\bar{Q}$ pair

The creation of a heavy quark, $Q$, involves a physical momentum transfer of
at least ${\cal O}(m_{Q})$, regardless of whether the production mechanism is
a hard process, MPI, or a $g\\!\rightarrow Q\bar{Q}$ shower branching. Given
an event-generator evolution algorithm that is ordered in a measure of
momentum transfer, when the evolution reaches a scale of order $m_{Q}$ one
could thus immediately veto any events that do not contain a specific desired
number of heavy quarks333Including not only quarks of flavour $Q$ but also
states that contain them, like $Q$-onia, or which can decay to them, e.g. via
$b\rightarrow c$ and/or $t\rightarrow b$ transitions.. Doing so effectively
“saves” the time it would otherwise have taken to evolve such events from
$m_{Q}$ to the perturbative cutoff, as well as the time required to hadronise
them. Since the running value of $\alpha_{s}$ is largest near the cutoff,
where the number of evolving partons is also highest, and hadronisation can
consume a further significant amount of time (in particular when colour
reconnections are involved), a veto at the scale $m_{Q}$ should produce order-
of-magnitude efficiency gains, compared with a standard simulation in which
events are fully evolved and hadronised before the event is inspected.

Figure 5: Differential creation probabilities for $b$ quarks, as a function of
the Pythia $p_{\perp}$ evolution parameter $p_{\perp\mathrm{evol}}$, for a
reference $gg\rightarrow gg$ hard process with
$\hat{p}_{\perp}=25~{}\mathrm{GeV}$ in proton-proton collisions at
$\sqrt{s}=10\,\mathrm{TeV}$. The solid red line shows FSR $g\rightarrow
b\bar{b}$ branchings, the dashed blue one shows ISR gluons backwards-evolving
to $b$ quarks, and the dotted black one shows MPI pair-creation and flavour-
excitation processes. The vertical dashed gray line indicates the default
value of the $b$ quark mass in Pythia, $m_{b}=4.8\,\mathrm{GeV}$.

A caveat to implementing such a strategy for Pythia, however, is that its
default transverse-momentum evolution scales [16, 17] for heavy-quark
processes are only guaranteed to be greater than $m_{Q}$ for hard processes
and for ISR, whereas the probability densities for FSR $g\\!\rightarrow
Q\bar{Q}$ splittings and MPI $gg\\!\rightarrow Q\bar{Q}$ processes are non-
zero all the way down to $p_{\perp\mathrm{evol}}\rightarrow 0$. This is
illustrated in Fig. 5, for a representative case of $gg\\!\rightarrow gg$
hard-scattering events with $\hat{p}_{\perp}=25\pm 0.5\,\mathrm{GeV}$ in
proton-proton collisions with $\sqrt{s}=10~{}\mathrm{TeV}$.

An interesting follow-up question for future work is thus whether it would be
physically justifiable to reformulate these algorithms in terms of a measure
that would associate all heavy-quark production processes with scales
$\geq~{}m_{Q}$. (This is, e.g., the choice made in the Vincia shower [21],
which however still relies on Pythia’s MPI evolution.)

For now, we accept that this subtlety will force a tradeoff between realising
the full possible efficiency gains and “missing” a small fraction of signal
events, which we will comment on further below. We therefore phrase our
implementation in terms of an arbitrary veto scale, not necessarily equal to
the heavy-quark mass, which can be varied to determine if an “acceptable”
tradeoff can be found, which may vary from application to application.
Specifically, we define the following Userhook vetoes:

Hard-process-level veto:

At this stage, nothing is known about the subsequent shower or MPI evolution,
except for what the starting scale for those evolutions will be. Our veto
function only accepts events that fulfil at least one of the following two
conditions: 1) the hard process itself contains the requisite heavy flavour
(by which we include any onium containing it or a heavier quark that can decay
to it), in which case a flag may also be set to bypass any downstream vetoes,
or 2) the starting scale for MPI and showers is above our user-defined veto
scale, so that we want to give MPI and/or showers a chance to produce the
heavy flavour. This essentially means that $gg\\!\rightarrow gg$ events with
$\hat{p}_{\perp}<{\cal O}(m_{Q})$ can be rejected already at this stage, with
minimum processing.

Event-evolution veto:

If the hard-scattering process did not contain the requisite heavy flavour but
was allowed a chance to produce it via MPI and/or showers, the event is
inspected again when the evolution reaches our veto scale, and is now rejected
if the required flavour (again including onia and/or heavier flavours) is
still not present in the event.

The improvement in efficiency when generating samples with these two UserHooks
is investigated for samples of events containing $b\bar{b}$ or $c\bar{c}$. The
time taken to generate the $Q\bar{Q}$ pairs is compared to a baseline without
the UserHooks included. All timing tests are performed using an Apple M1
MacBook Pro.444The timing studies were performed using single-core jobs.
Benchmarking tests suggest in this configuration the machine has a CPU power
of approximately 44 HS06. The relative speed-up and fraction of events missed
due to the evolution scale definition are shown for $b\bar{b}$ pairs in Fig.
6.

Figure 6: (Left) Relative speed enhancements of Pythia when generating
$b\bar{b}$ events with process-level and evolution-level UserHooks. (Right)
The fractions of $b\bar{b}$ events that are not retained by the UserHooks as a
function of $\hat{p}_{T}$ scale.

A significant improvement in efficiency is found when generating $b\bar{b}$
pairs with the UserHooks. The improvement is less significant when generating
$c\bar{c}$ pairs because the smaller $c$-quark mass means the event evolution
must continue further before the event can be vetoed.

The $p_{T}$ distribution of $B$ hadrons in events that are not retained by the
UserHooks are shown in Fig. 7. This sample, produced with the Simple Shower
model misses ${b}{\overline{b}}$ pairs produced in both the parton shower and
as additional MPI interactions. Overall, setting a $\hat{p}_{T}$ scale of
$4\mathrm{\,Ge\kern-1.00006ptV}$ gives a factor 10 improvement in simulation
speed, but leads to a small distortion in the $p_{T}$ spectra of the generated
$b$ hadrons.

Figure 7: Kinematic distribution of $b$-hadrons in Pythia samples generated
with the Simple Shower model. Those that are accepted or missed by the
UserHooks discussed in the text with a $\hat{p}_{T}$ scale of
$4\mathrm{\,Ge\kern-0.90005ptV}$ are highlighted and further split according
to whether the missed heavy quarks originated during the parton shower or as
an additional MPI process.

The impact of these efficiency improvements can further contextualised in
terms of the typical time taken to generate specific singly-heavy hadrons with
Pythia. The typical times with and without the developed UserHooks are listed
in Table 1.

Hadron | No UH | With UH
---|---|---
$B^{+}$ | 0.15 | 0.013
$B^{0}$ | 0.15 | 0.013
$B_{s}^{0}$ | 0.55 | 0.045
$\Lambda_{b}^{0}$ | 1.5 | 0.13
$\Sigma_{b}^{-}$ | 27 | 1.9
$\Sigma_{b}^{0}$ | 26 | 2.3
$\Sigma_{b}^{+}$ | 26 | 2.3
$\Xi_{b}^{-}$ | 11 | 0.9
$\Xi_{b}^{0}$ | 13 | 0.8
${\Xi^{\prime}}_{b}^{-}$ | 220 | 18
${\Xi^{\prime}}_{b}^{0}$ | 150 | 17
$\Omega_{b}^{-}$ | 310 | 28

Hadron | No UH | With UH
---|---|---
$D^{0}$ | 0.011 | 0.006
$D^{+}$ | 0.017 | 0.009
$D_{s}^{+}$ | 0.04 | 0.02
$\Lambda_{c}^{0}$ | 0.11 | 0.06
$\Sigma_{c}^{++}$ | 2.2 | 1.2
$\Sigma_{c}^{+}$ | 2.0 | 1.0
$\Xi_{c}^{+}$ | 0.7 | 0.4
$\Xi_{c}^{0}$ | 0.8 | 0.4
${\Xi^{\prime}}_{c}^{+}$ | 11 | 6
${\Xi^{\prime}}_{c}^{0}$ | 12 | 9
$\Omega_{c}^{0}$ | 36 | 12

Table 1: Typical generation times in seconds for single-heavy hadrons with and
without the Userhooks (UH) described in the text. The tests are performed
using the Simple Shower model. The vetoes are imposed with $\hat{p}_{T}$ cut-
off scales of $4.0{\mathrm{\,Ge\kern-0.90005ptV\\!/}c}$ and
$1.5{\mathrm{\,Ge\kern-0.90005ptV\\!/}c}$ for the $b\bar{b}$ and $c\bar{c}$
samples, respectively. Due to the size of the generated samples the
uncertainty in the typical generation time for $\Sigma$, $\Xi$ and $\Omega$
baryons is up to 7%, 15% and 30%, respectively. The corresponding
uncertainties for the meson and $\Lambda$ samples are negligible.

###  4.2 Simulating final states involving multiple $Q\bar{Q}$ pairs

When simulating events with multiple pairs of heavy quarks, the same Userhook
vetoes described previously in Section 4.1 can be utilised. In principle, for
$B_{c}^{+}$ production and the like, it would be useful to apply an event-
evolution veto first at ${\cal O}(m_{b})$ and then again at ${\cal O}(m_{c})$.
However, with current versions of Pythia, the event can only be inspected at a
single value of the evolution scale. Therefore, we set the the event-evolution
threshold according to the heaviest quark being simulated, while the end-of-
evolution parton-level veto is used to check for any required secondary heavy
flavour.

Parton-level veto:

If both a $b\bar{b}$ and a $c\bar{c}$ pair is requested (and/or onia
containing them), the presence of the lighter of the two flavours is checked
for at the end of the parton-level evolution.

We believe that, by making minor modifications to the Pythia source code it
could be possible to allow the event-evolution veto to inspect the event
multiple times, at different scales during the event evolution (e.g. at both
the $b$ and $c$ quark masses). This would save the time spent evolving the
event from the lower mass scale to the hadronisation scale, providing benefits
in addition to those demonstrated here.

###  4.3 Simulating final states with specific hadrons

Further UserHooks can be placed at the partonic event level to increase the
efficiency of generating specific hadrons. For example, when generating
$B_{c}^{+}$ mesons, the invariant mass of all $c$ and $\bar{b}$ pairs (and
vice versa) can be calculated to veto events in which there are no
combinations that have a small enough invariant mass to form a hadron. The
timing improvements when using the additional levels of UserHook are shown in
Fig. 8, where the first UserHook to be included requires the smallest
invariant mass of any $b\bar{c}$ or $\bar{b}c$ pair to be below
$10{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, the second requires the
partonic-level event to contain both $b\bar{b}$ and $c\bar{c}$ pairs, and the
third and fourth correspond to the hard process and event evolution vetoes
already described, here applied just to the $b$ quark.

Figure 8: Relative speed enhancements of Pythia when generating $B_{c}^{+}$
mesons with the different UserHooks included. Each point includes an
additional UserHook as described in the text.

Using combinations of the UserHooks described above it is possible to improve
the efficiency of generating both singly-heavy and doubly-heavy hadrons with
Pythia. Typical times to generate doubly-heavy hadrons are listed in Table 2.

The UserHooks that inspect the event at the parton-level or after do not
result in any further missed hadrons, in contrast to the hard-process level
and event-evolution vetoes discussed in Section 4.1. Therefore when generating
doubly-heavy hadron with these extra requirements, the fraction and kinematic
distributions of the missed hadrons will be the same as for singly-heavy
hadrons.

Hadron |  UserHook requirement | No UH | With UH | Speed-up
---|---|---|---|---
${B}_{c}^{+}$ |  ${b}{\overline{b}}{c}{\overline{c}}$ and $m({b}{\overline{c}})<10{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ | 47 | 3.4 | $\sim 14$
${\mathchar 28932\relax}^{++}_{{c}{c}}$ | ${c}{\overline{c}}{c}{\overline{c}}$ | 84 | 37 |
${\mathchar 28932\relax}^{+}_{{c}{c}}$ | ${c}{\overline{c}}{c}{\overline{c}}$ | 78 | 37 | $\sim 2$
$\Omega_{cc}^{0}$ | ${c}{\overline{c}}{c}{\overline{c}}$ | 620 | 260 |
$\Xi_{bc}^{+}$ | ${b}{\overline{b}}{c}{\overline{c}}$ | 830 | 36 |
$\Xi_{bc}^{0}$ | ${b}{\overline{b}}{c}{\overline{c}}$ | 900 | 41 | $\sim 20$
$\Omega_{bc}^{0}$ | ${b}{\overline{b}}{c}{\overline{c}}$ | 3300 | 300 |
Table 2: Typical generation times in seconds and the approximate speed-up
factors for doubly-heavy hadrons with and without UserHooks. The minimal quark
requirements used for each UserHook are included, along with any additional
requirements. Note, the baryonic samples are generated with a different colour
reconnection model described later in Section 5.2. For comparison, generating
${B}_{c}^{+}$ mesons using Bcvegpy takes approximately 4${\mathrm{\,ms}}$
seconds to generate the hard process and a further 4${\mathrm{\,ms}}$ to
perform the full event evolution with Pythia.

##  5 Comparison of doubly-heavy hadron kinematics

The efficiency improvements detailed in the previous section enable a range of
comparisons that were previously unfeasible. For example, 250 000
${B}_{c}^{+}$ mesons have been generated for this study which would take
approximately 10 CPU days with the UserHooks, or 140 CPU days without. In this
section distributions of doubly-heavy hadrons generated with Pythia are
compared to other standalone generators.

###  5.1 $B_{c}^{+}$ mesons

The production of ${B}_{c}^{+}$ mesons at hadron colliders is generally
calculated in terms of the SPS process
$gg\\!\rightarrow{{B}_{c}^{+}}{b}{\overline{c}}$ that is assumed to dominate
[22, 23, 24], although there are some suggestions that flavour excitation
processes could play an important role [25]. Currently, DPS production has not
been considered. Measurements of the ratio of the ${B}_{c}^{+}$ to the
${{B}^{+}}$ cross-section times branching fraction are found to be consistent
with predictions of SPS contributions only [26], although the theoretical
predictions have large uncertainties. In contrast, measurements of $\Upsilon
D$ cross section are consistent with the presence of DPS [4]. The $\Upsilon D$
system also requires the creation of a ${b}{\overline{b}}$ and a
${c}{\overline{c}}$ pair, implying that there is the potential for
${B}_{c}^{+}$ mesons to be similarly produced.

Samples of ${B}_{c}^{+}$ mesons are generated with Pythia and Bcvegpy. The
specific generation settings are listed in Table 3.

Pythia option | Setting
---|---
Beams:eCM | 13000
SoftQCD:nonDiffractive | on
PartonShowers:Model | 1
PartonLevel:MPI | on

Bcvegpy option | Setting
---|---
ENERGYOFLHC | 13000
naccel | 2
i_mix | 1
mix_type | 1
igenmode | 1
istate | 1
numofevents | 100000
ncall | 10000
n_itmx | 15
$m({b})$ | 5$\mathrm{\,Ge\kern-1.00006ptV}$
$m({c})$ | 1.275$\mathrm{\,Ge\kern-1.00006ptV}$

Table 3: Pythia and Bcvegpy generation settings used to simulate the samples
of ${B}_{c}^{+}$ mesons. The specific meanings of the parameters can be found
in the relevant documentation [12, 13]. Additionally, in Pythia all ground-
state heavy hadrons are treated as stable. As such, contributions from $X_{c}$
hadrons originating from $X_{b}$ decays are explicitly excluded from the
samples.

The kinematic distributions and differential cross-sections are compared in
Fig. 9. The sample generated with Pythia is split according to the origin of
the $c$-quark and $\overline{b}$-quark that formed the ${B}_{c}^{+}$ meson:
those that originated in the same parton-parton interaction are categorised as
SPS and those from different parton-parton interactions are categorised as
DPS. This is determined by comparing the parent quarks to the record of parton
interactions stored in Pythia’s PartonSystems class. The production cross-
section predicted by Pythia receives a significant contribution from DPS and
is found to be larger than the production cross-section obtained from Bcvegpy.

Figure 9: Kinematic distributions of $B_{c}^{+}$ mesons generated with Bcvegpy
and Pythia.

The kinematic distributions of the ${B}_{c}^{+}$ alone do not provide
significant discrimination power to the presence of hadron formation in DPS,
other than in the absolute normalisation of the cross section. Measurements of
the ${B}_{c}^{+}$ cross section relative to ${{B}^{+}}$ mesons have been
performed [26, 27, 28], but extracting the absolute cross section requires
theoretical predictions of the branching fractions.

A better way to investigate the role of SPS and DPS contributions is to
measure the production cross section as a function of the number of parton-
parton interactions in a collision. For hadrons formed in SPS processes,
increasing the number of parton interactions would linearly increase the
number of opportunities to form the hadron, as each new parton interaction
would present one more opportunity for the hadron to form. However, hadrons
formed in DPS processes would see the rate of formation increase quadratically
with the number of interactions, as each hadron requires two parton
interactions to form. These different relationships can be exploited to
differentiate the components by considering the ratio of doubly-heavy to
singly-heavy hadron cross sections, as a function of the number of parton-
parton interactions. This ratio would be flat if singly- and doubly-heavy
hadrons are produced by the same mechanism — SPS — while it would increase
linearly if there is a nontrivial DPS component to doubly-heavy hadron
production.

Figure 10: Ratio of differential cross-sections of $B_{c}^{+}$ and ${{B}^{+}}$
mesons as a function of (top left) the number of parton-parton interactions in
a collision and (top right) the number of charged particles within the pseudo-
rapidity region $2.0<\eta<4.5$, as generated with Bcvegpy and Pythia.
Uncertainties are from simulation statistics only.

In Pythia, both mechanisms are present, while in Bcvegpy, a single
$gg\rightarrow{{B}_{c}^{+}}{b}{\overline{c}}$ interaction is produced for each
event which is then passed to Pythia for showering, MPI, and hadronisation. In
this case, there is therefore no opportunity for heavy quarks from different
parton-parton interactions to form the ${B}_{c}^{+}$ meson and the production
is independent of the total number of parton-parton interactions.

The cross-section ratio of ${B}_{c}^{+}$ to ${{B}^{+}}$ mesons is compared for
Pythia and Bcvegpy in Fig. 10 as a function of the number of parton-parton
interactions. In this figure no kinematic requirements have been placed on the
rapidity or transverse momentum of the ${B}_{c}^{+}$ meson or final-state
particles. As expected, the contribution from DPS varies as a function of the
number of parton interactions in the event. A significant enhancement is seen
in events with many parton interactions. This quantity is not directly
observable, but it is highly correlated to the number of particles in a
collision. A more realistic, experimentally-accessible distribution is
therefore the cross-section ratio as a function of the number of charged
particles within a given acceptance. In the right pane of Fig. 10 and
throughout we use the acceptance of the LHCb experiment ($2.0<\eta<4.5$), but
other ranges would show similar results.

Figure 11:  Distribution of the number of charged particles within
$2.0<\eta<4.5$ in events containing a ${B}_{c}^{+}$ meson as generated with
Bcvegpy and Pythia.

The distribution of the number of charged particles within $2.0<\eta<4.5$ in
events containing a ${B}_{c}^{+}$ meson is shown in Fig. 11 for the different
generation configurations. In the Pythia SPS and Bcvegpy samples there are
means of approximately 34 charged particles per event, whilst the Pythia DPS
sample has a higher mean of around 41.

The significant difference in these distributions would allow measurements of
the differential cross section with respect to ${{B}^{+}}$ mesons to
differentiate between ${B}_{c}^{+}$ mesons produced in SPS and DPS.

###  5.2 $\Xi_{cc}^{++}$ baryons

The doubly-charmed baryon ${\mathchar 28932\relax}^{++}_{{c}{c}}$ has been
observed [29, 30], and could similarly receive contributions from DPS
production mechanisms. Samples of $\Xi_{cc}^{++}$ baryons are simulated with
both Pythia and GenXicc. In Pythia, the formation of baryons with multiple
heavy quarks is only possible with the so-called QCD colour-reconnection
option [8] which allows for the formation of doubly-heavy diquarks via string
junctions [31]. The specific Pythia settings are listed in Table 4.

Pythia option | Setting
---|---
Beams:eCM | 13000
SoftQCD:nonDiffractive | on
PartonShowers:Model | 1
PartonLevel:MPI | on
StringPT:sigma | 0.335
StringZ:aLund | 0.36
StringZ:bLund | 0.56
StringFlav:probQQtoQ | 0.078
StringFlav:ProbStoUD | 0.2
StringFlav:probQQ1toQQ0join | 0.0275, 0.0275, 0.0275, 0.0275
MultiPartonInteractions:pT0Ref | 2.15
BeamRemnants:remnantMode | 1
BeamRemnants:saturation | 5
ColourReconnection:mode | 1
ColourReconnection:allowDoubleJunRem | off
ColourReconnection:m0 | 0.3
ColourReconnection:allowJunctions | on
ColourReconnection:junctionCorrection | 1.20
ColourReconnection:timeDilationMode | 2
ColourReconnection:timeDilationPar | 0.18
Table 4: Pythia generation settings used to simulate the samples of
${\mathchar 28932\relax}^{++}_{{c}{c}}$ mesons.

The GenXicc generator can produce different configurations of the initial
charm hadrons. The processes $gg\\!\rightarrow{{\mathchar
28932\relax}^{++}_{{c}{c}}}{\overline{c}}{\overline{c}}$ and
$gc\\!\rightarrow{{\mathchar 28932\relax}^{++}_{{c}{c}}}{\overline{c}}$ are
simulated. The kinematic distributions of the ${\mathchar
28932\relax}^{++}_{{c}{c}}$ samples are shown in Fig 12, where again the
Pythia samples have been split according to the whether the two heavy quarks
originated from the same or different parton-parton interactions. Analogously
to the case for ${B}_{c}^{+}$, a significantly larger cross section is
predicted by Pythia, due to the presence of the DPS component.

Figure 12: Kinematic distributions of $\Xi_{cc}^{++}$ baryons generated with
Pythia and GenXicc.

Figure 13: Ratio of differential cross-section of ${\mathchar
28932\relax}^{++}_{{c}{c}}$ and ${D}^{+}$ hadrons as a function of (left) the
number of parton interactions in a collision and (right) the number of charged
particles within the pseudo-rapidity region $2.0<\eta<4.5$, as generated with
GenXicc and Pythia.

The ratio of the differential cross-section between the ${\mathchar
28932\relax}^{++}_{{c}{c}}$ and ${D}^{+}$ hadrons are shown as a function of
the number of parton interactions and number of charged particles within the
pseudo-rapidity region $2.0<\eta<4.5$ in Fig. 13. The cross-section ratio is
clearly able to differentiate between DPS and SPS production mechanisms. The
difference in the cross-section ratios for the $gg\\!\rightarrow{{\mathchar
28932\relax}^{++}_{{c}{c}}}{\overline{c}}{\overline{c}}$ and
$gc\\!\rightarrow{{\mathchar 28932\relax}^{++}_{{c}{c}}}{\overline{c}}$
processes simulated by GenXicc may be caused by the different typical momentum
transfers in the two processes. When simulating the rest of the underlying
event with Pythia, this quantity determines the maximum scale of the event
evolution, therefore the typically lower momentum transfers in
$gc\\!\rightarrow{{\mathchar 28932\relax}^{++}_{{c}{c}}}{\overline{c}}$
processes may cause fewer parton-parton interactions to be generated in the
subsequent event evolution.

###  5.3 $\Xi_{bc}^{+}$ baryons

Similarly, samples of $\Xi_{bc}^{+}$ can be produced with Pythia and GenXicc.
These baryons have not yet been observed. The kinematic distributions of the
$\Xi_{bc}^{+}$ baryons in samples generated with Pythia and GenXicc are
compared in Fig 14, and the ratios of the differential cross section with
respect to the ${{B}^{+}}$ meson are shown in Fig. 15 as a function of both
the number of parton interactions and number of charged tracks within
$2.0<\eta<4.5$. The production in Pythia is again dominated by a significant
DPS contribution.

Figure 14: Kinematic distributions of $\Xi_{bc}^{+}$ baryons generated with
Pythia and GenXicc.

Figure 15: Differential cross-section of $\Xi_{bc}^{+}$ baryons as a function
of (left) the number of parton interactions in a collision and (right) the
number of charged tracks within the rapidity region $2.0<y<4.5$, as generated
with GenXicc and Pythia.

###  5.4 Quarkonia

For the purposes of this study, quarkonia are considered doubly-heavy hadrons
as they contain two heavy quarks. However, unlike the previously discussed
doubly-heavy hadrons, quarkonia can be formed from a single heavy-quark pair.
As a result, contributions to quarkonium production from heavy quarks in
different parton interactions will only give a subleading contribution to the
total production rate. This can be observed in Fig. 16 where the ratio of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ to ${D}^{+}$ production cross
sections are shown as a function of the number of parton interactions and the
number of charged tracks within $2.0<\eta<4.5$. There is a contribution from
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mesons formed from $c$ and
$\overline{c}$ quarks from different parton interactions, however this is much
smaller than the total rate of ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
production. To increase the fraction of events with DPS contributions, events
with both a ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson and two
additional charm hadrons can be reconstructed. This removes events in which
there is only a single quark pair and leads to a different set of measurements
that can be made as discussed in the next Section.

Figure 16: Ratio of differential cross-section of
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ and ${D}^{+}$ hadrons as a
function of (left) the number of parton interactions in a collision and
(right) the number of charged particles within the pseudo-rapidity region
$2.0<\eta<4.5$, as generated with Pythia.

##  6 Associated production of singly- and doubly-heavy hadrons in events with
multiple $Q\bar{Q}$ pairs

In proton-proton collisions that produce two pairs of heavy quarks, i.e.
$c\bar{c}c\bar{c}$, $c\bar{c}b\bar{b}$ or $b\bar{b}b\bar{b}$, information
about the production mechanisms can be inferred from the relative properties
of a doubly-heavy hadron and two singly-heavy hadrons that can be formed from
the additional heavy quarks. Examples of the different combinations of doubly-
and singly-heavy associated productions are listed in Table 5, where singly-
heavy hadrons containing a heavy quark are represented by $X_{Q}$. Only the
${B}_{c}^{+}$ meson, $b$\- and $c$-flavoured quarkonia and ${\mathchar
28932\relax}^{++}_{{c}{c}}$ baryon have currently been observed [32, 33, 34,
35, 29, 30].

| $c\bar{c}c\bar{c}$ | $c\bar{c}b\bar{b}$ | $b\bar{b}b\bar{b}$
---|---|---|---
Doubly-heavy meson |  | $B_{c}^{+}X_{b}X_{\bar{c}}$ |
(excluding quarkonia) |  |  |
Doubly-heavy baryons | $\Xi_{cc}X_{\bar{c}}X_{\bar{c}}$ | $\Xi_{bc}X_{\bar{b}}X_{\bar{c}}$ | $\Xi_{bb}X_{\bar{b}}X_{\bar{b}}$
| $\Omega_{cc}X_{\bar{c}}X_{\bar{c}}$ | $\Omega_{bc}X_{\bar{b}}X_{\bar{c}}$ | $\Omega_{bb}X_{\bar{b}}X_{\bar{b}}$
Quarkonia | $\psi(nS)X_{c}X_{\bar{c}}$ | $\psi(nS)X_{b}X_{\bar{b}}$ | $\Upsilon(nS)X_{b}X_{\bar{b}}$
|  | $\Upsilon(nS)X_{c}X_{\bar{c}}$ |
Table 5: Examples of combinations of doubly- and singly-heavy hadrons in
processes with two pairs of heavy quarks.

Of the processes listed in Table 5, those categorised as doubly-heavy mesons
and doubly-heavy baryons require quarks from different $Q\bar{Q}$ pairs to
combine in order to form. This configuration is referred to as mixed. In
contrast, the quarkonium in $\psi(nS)X_{b}X_{\bar{b}}$ or
$\Upsilon(nS)X_{c}X_{\bar{c}}$ events contain differently flavoured quarks
compared to the associated singly-heavy hadrons, and therefore must be formed
from an individual heavy-quark pair. This configuration is referred to as
unmixed. In events containing a quarkonium, $X_{Q}$ and $X_{\bar{Q}}$ hadrons
where the flavours of the singly-heavy hadrons match the doubly-heavy flavour,
the quarkonium could have been formed from a single $Q\bar{Q}$ pair, i.e.
unmixed, or from the combination of a $Q$ and $\bar{Q}$ from two different
$Q\bar{Q}$ pairs, i.e. mixed. This is demonstrated for the charm quark in Fig.
17.

(a) Mixed

(b) Unmixed

Figure 17: Possible configurations of quark lines in DPS events with two pairs
of charm quarks.

Therefore, for the ${B}_{c}^{+}$, $\Xi_{cc}$ and $\Xi_{bc}$ hadrons any
indication that DPS processes contribute to the production would indicate that
heavy quarks from different parton interactions can form hadrons. In contrast,
the indications that DPS processes contribution to $\psi(nS)X_{b}X_{\bar{b}}$
or $\Upsilon(nS)X_{c}X_{\bar{c}}$ events implies that there are events with
multiple heavy-quark pairs from different parton interactions, but doesn’t
give any information about the hadronisation of those different pairs. The
quarkonia final states $\psi(nS)X_{c}X_{\bar{c}}$ and
$\Upsilon(nS)X_{b}X_{\bar{b}}$ can give information about both the presense of
multiple heavy quark pairs and hadronisation in DPS, but care must be taken to
separate the mixed and unmixed contributions.

The kinematic distributions of the singly- and doubly-heavy hadrons depend on
the specific production mechanism and differ for those possible within SPS or
DPS. Collisions in which both heavy hadron pairs originate from a SPS, for
example the process shown in Fig. 3(a), have correlations between the
kinematics of at least one pair of heavy quarks. Processes that involve DPS
introduce the possibility of double flavour excitation processes, such as Fig.
3(c), in which the kinematics of the heavy quarks that form the hadron can be
relatively uncorrelated to the remaining two companion heavy quarks.555The
small fraction of hadrons that are missed by the UserHooks (Sec. 4.1) are not
found to significantly alter the relative kinematic correlations.

The production of quarkonia is complicated further by contributions from
colour-octet mechanisms [36, 37, 38], which are included by default in the
Pythia simulation samples.

###  6.1 Studies with $B_{c}^{+}X_{b}X_{\bar{c}}$ events

The relative kinematic distribution in $B_{c}^{+}X_{b}X_{\bar{c}}$ events are
studied to determine if there are differences between the DPS and SPS sub-
samples for the mixed configuration, and to compare with the standalone
generator Bcvegpy. The angular separation in the transverse view between the
${B}_{c}^{+}$ and the associated $X_{b}$ hadron
$\Delta\phi({{B}_{c}^{+}},X_{b})$ is plotted against the same quantity between
the ${B}_{c}^{+}$ and associated $X_{{\overline{c}}}$ hadron
$\Delta\phi({{B}_{c}^{+}},X_{\bar{c}})$ in Fig. 18 for the different
generators and configurations. To ensure an unambiguous association between
the heavy quarks, only events with a total of one ${B}_{c}^{+}$ meson, one
$X_{\bar{c}}$ hadron and one $X_{b}$ hadron are used.

(a) Bcvegpy

(b) Pythia\- Just DPS

(c) Pythia\- Just SPS

Figure 18: Angular separation in the transverse view between $B_{c}^{+}$
mesons and the associated $X_{b}$ or $X_{\bar{c}}$ hadron in events generated
with the Bcvegpy and Pythia generators.

In the sample generated with Bcvegpy the generated events are found
predominately with large angles between the $B_{c}^{+}$ meson and $X_{b}$
hadron, and small angles between the $B_{c}^{+}$ and $X_{\bar{c}}$ hadrons.
This is consistent with the topology shown in Fig. 3(a) in which the
$b\bar{b}$ pair are produced in the hardest interaction, and therefore back-
to-back in the transverse plane. The $c$-quark resulting from a gluon
splitting would then be produced in a cone around the $B_{c}^{+}$ direction.
In contrast, the DPS sample produced by Pythia introduces the possibility of
additional production mechanisms including those shown in Fig. 3(b) and Fig.
3(c). As such the distribution of events in the 2D plane is less localised as
a result of the contributions from many different associated production
correlations.

###  6.2 Studies with ${\mathchar 28935\relax{(1S)}}X_{c}X_{\bar{c}}$ events

Events containing both an ${\mathchar 28935\relax{(1S)}}$ meson and
$X_{c}X_{\bar{c}}$ pair can only receive contributions from the configuration
referred to as unmixed. It therefore provides a suitable system to test
whether MPIs contribute significantly to events with multiple pairs of heavy
quarks, but cannot provide insight into the hadronisation of heavy quarks from
different parton interactions.

(a) Pythia\- Just SPS

(b) Pythia\- Just DPS

(c)

(d)

Figure 19: (Top) Relative transverse distributions in ${\mathchar
28935\relax{(1S)}}X_{c}X_{\bar{c}}$ events for SPS and DPS processes. (Bottom
left) One-dimensional projections of the relative angular distributions,
compared to measurements from Ref. [4]. (Bottom right) Relative differential
cross-section of ${\mathchar 28935\relax{(1S)}}X_{c}X_{\bar{c}}$ events with
respect to ${{B}^{+}}$ as a function of the number of charged particles within
the pseudo-rapidity region $2.0<\eta<4.5$.

The relative transverse distributions between the ${\mathchar
28935\relax{(1S)}}$ meson and the $X_{c}$ and $X_{\bar{c}}$ hadrons are shown
in Fig. 19(a) and 19(b) . The two-dimensional distributions show a clear
difference in the relative distributions of the two types of process. In the
SPS process the ${\mathchar 28935\relax{(1S)}}$ meson has a strong tendency to
be produced back-to-back to both charm hadrons in the transverse plane. This
could result from $gg\rightarrow{\mathchar 28935\relax{(1S)}}g$ parton
interactions in which the outgoing gluon subsequently splits
$g\rightarrow{c}{\overline{c}}$. In contrast, in DPS processes there is little
correlation between the transverse directions. The predictions can be compared
to measurements by LHCb [4] by making a one-dimensional projection, as shown
in Fig. 19(c). Only events in which both the ${\mathchar 28935\relax{(1S)}}$
and one of the corresponding $X_{c}$ hadron are within the LHCb acceptance are
compared in this figure. The data are consistent with the predictions
including DPS, but the significant difference between the SPS and DPS samples
is diluted when projected onto one dimension. This strongly motivates
performing measurements in which both associated hadrons are reconstructed,
such that the two dimensional distributions can be determined. Finally, the
ratio of the differential cross section relative to the ${{B}^{+}}$ meson is
shown in Fig. 19(d) as a function of the number of charged particles within
$2.0<\eta<4.5$. It similarly shows a strong separation power between the two
process types.

###  6.3 Studies with $J/\psi X_{c}X_{\bar{c}}$ events

Samples containing a ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson and
two charm hadrons are generated with Pythia. The samples are split according
to whether the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson was
created in a single interaction, or from quarks in different parton
interactions, as before. Additionally, the ancestors of the charm quarks are
studied to determine if the two quark lines are mixed or unmixed, as defined
in Fig. 17.

Figure 20: Ratio of differential cross-section of
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X_{{c}}X_{{\overline{c}}}$ and
${D}^{+}$ hadrons as a function of (left) the number of parton interactions in
a collision and (right) the number of charged particles within the pseudo-
rapidity region $2.0<\eta<4.5$, as generated with Pythia.

The ratio of the differential cross sections with respect to the ${D}^{+}$
meson is compared for the four categories, as shown in Fig. 20. The two DPS
categories show an increasing trend with both the number of parton
interactions and charged particles within $2.0<\eta<4.5$. The fraction of
events produced in DPS processes is significantly increased with respect to
the samples in which just the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
was reconstructed in Section 5.4.

The relative transverse directions of the three particles is shown for the
four categories in Fig. 21. The two unmixed categories show similar features
to those already seen in ${\mathchar
28935\relax{(1S)}}X_{{c}}X_{{\overline{c}}}$ events. In contrast the two mixed
categories have distinct distributions, with the DPS mixed category having a
slight tendency to have the two charm quarks produced in the same direction as
the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$, a topology that may arise
due to events with two flavour excitation interactions.

(a) Pythia\- SPS unmixed

(b) Pythia\- DPS unmixed

(c) Pythia\- SPS mixed

(d) Pythia\- DPS mixed

Figure 21: Relative transverse distributions in
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X_{c}X_{\bar{c}}$ events for
SPS and DPS processes in both the mixed and unmixed configurations.

The presence of these events categorised as DPS mixed would indicate that
heavy quarks from different parton interactions can hadronise to form a
$J/\psi$ meson. Using a combination information from the relative cross
section ratio and relative transverse directions, the presence of this
contribution could be determined.

###  6.4 Studies with ${\mathchar 28935\relax{(1S)}}X_{b}X_{\bar{b}}$ events

Similar studies can be performed on the ${\mathchar
28935\relax{(1S)}}X_{{b}}X_{{\overline{b}}}$ system. Large samples are
generated with Pythia requiring a minimal event content of
${b}{\overline{b}}{b}{\overline{b}}$. Events with one ${\mathchar
28935\relax{(1S)}}$ meson and a $X_{b}X_{\overline{b}}$ pair are selected, and
the relative distributions studied. The events are split according to whether
the $b$ and $\overline{b}$ that formed the ${\mathchar 28935\relax{(1S)}}$
were from the same or different parton systems. As before, the ancestry of the
$b$ and $\overline{b}$ quarks is studied to determine if the configuration is
mixed or unmixed. The relative differential cross-section relative to the
${{B}^{+}}$ meson is shown in Fig. 22 as a function of the number of parton
interactions and the number of charged particles within $2.0<\eta<4.5$. A
similar trend is observed in this system.

Figure 22: Ratio of differential cross-section of ${\mathchar
28935\relax{(1S)}}X_{{b}}X_{{\overline{b}}}$ and ${{B}^{+}}$ hadrons as a
function of (left) the number of parton interactions in a collision and
(right) the number of charged particles within the pseudo-rapidity region
$2.0<\eta<4.5$, as generated with Pythia.

Additionally, the relative transverse directions of the ${\mathchar
28935\relax{(1S)}}$, $X_{{b}}$ and $X_{{\overline{b}}}$ hadrons is studied for
the four production categories. The distributions, shown in Fig. 23 show a
similar distribution for the SPS unmixed as before. In this system, events
categorised as DPS unmixed have a slight tendency to have an anti-correlation
between $\Delta\phi({\mathchar 28935\relax{(1S)}},X_{b})$ and
$\Delta\phi({\mathchar 28935\relax{(1S)}},X_{\bar{b}})$. This could correspond
to events in which one parton interaction produces the ${\mathchar
28935\relax{(1S)}}$ meson and another produces a $b$ $\overline{b}$ pair from
pair production that are back-to-back in the transverse plane. Therefore as
$\Delta\phi({\mathchar 28935\relax{(1S)}},X_{b})$ increases, the corresponding
value of $\Delta\phi({\mathchar 28935\relax{(1S)}},X_{\bar{b}})$ decreases.
The DPS unmixed distribution has a distinct distribution, potentially allowing
discrimination.

(a) Pythia\- SPS unmixed

(b) Pythia\- DPS unmixed

(c) Pythia\- SPS mixed

(d) Pythia\- DPS mixed

Figure 23: Relative transverse distributions in ${\mathchar
28935\relax{(1S)}}X_{b}X_{\bar{b}}$ events for SPS and DPS processes in both
the mixed and unmixed configurations.

###  6.5 Impact of the Colour Reconnection model

As mentioned briefly in Sec. 5.2 and discussed in more detail in [8, 9, 12],
the QCD CR model of [8] allows for coherent combinations of three SU(3)
triplets (or antitriplets) into colour-singlet states, connected via Y-shaped
“string junctions” [31]. In the limit of a small invariant mass between two
heavy quarks, the junction topology (when allowed by the CR selection rules
[8]) reduces to a doubly-heavy diquark. In Pythia, this mechanism is
essentially the sole means by which doubly-heavy baryons can be produced at
all, making the rates and spectra of such baryons — and their dependence on
event characteristics such as $N_{\mathrm{Charged}}$ — particularly sensitive
probes of this type of colour-space ambiguities. (Accordingly, we note that
all of our Pythia results for doubly-heavy baryons in this work were produced
with the QCD CR option; otherwise the rates would be essentially zero.)

The choice of CR model also has an impact on the relative size of the DPS
contribution to doubly-heavy meson production. The effect of using the QCD CR
model described in Table 4 on the slope of the ratio of ${B}_{c}^{+}$ to
${{B}^{+}}$ differential cross sections is shown in Fig 24. With the QCD CR
model, the DPS contribution gets smaller, altering the slope of the total
distribution. As such, measurements of this distribution may help
differentiate between the different models.

Figure 24: Ratio of differential cross-section of ${B}_{c}^{+}$ and
${{B}^{+}}$ hadrons as a function of (left) the number of parton interactions
in a collision and (right) the number of charged particles within the pseudo-
rapidity region $2.0<\eta<4.5$, as generated with Pythia both with the
alternative CR options specified in Table 4 (triangles) and with the default
CR options (circles).

Similarly, the relative cross section distributions are found to vary in
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X_{c}X_{\bar{c}}$ samples
generated with the QCD CR options enabled, as shown in Fig. 25. The
contribution from the DPS mixed configuration decreases whilst DPS unmixed
increases, implying that the likelihood for the heavy quarks from different
parton-parton interactions to be combined into a single hadron is sensitive to
the choice of CR scheme.

Figure 25: Ratio of differential cross-section of
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X_{{c}}X_{{\overline{c}}}$ and
${D}^{+}$ hadrons as a function of (left) the number of parton interactions in
a collision and (right) the number of charged particles within the pseudo-
rapidity region $2.0<\eta<4.5$, as generated with Pythia both with the
alternative CR options specified in Table 4 (triangles) and with the default
CR options (circles).

Additionally, the relative transverse direction distributions are found to
differ for the DPS configurations. The corresponding distributions are shown
in Fig. 26. Measurements of $J\/\psi X_{c}X_{\bar{c}}$ events may additionally
help to differentiate between the different models of colour reconnection.

(a) Pythia\- DPS unmixed

(b) Pythia\- DPS mixed

Figure 26: Relative transverse distributions in
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X_{c}X_{\bar{c}}$ events for
DPS processes in both the mixed and unmixed configurations, as generated with
Pythia with the colour reconnection options specified in Table 4.

##  7 Experimental measurements and feasibility

In the Pythia simulation studies performed for this paper, the production of
doubly-heavy hadrons is predicted to have a significant contribution from DPS
production processes. New measurements of the relative cross section for the
doubly-heavy hadrons with respect to singly-heavy hadron as a function of the
collision multiplicity would help identify if such contributions are present
in nature, as proposed in Section 5. Unlike recent observations of strangeness
enhancements in the ratio of ${B}^{0}_{s}$ to ${B}^{0}$ cross sections [39],
the enhancements from DPS are not expected to be localised. The most suitable
doubly-heavy hadron for this would be the ${B}_{c}^{+}$ meson. The significant
yields reported in a selection of different papers are listed in Table 6.
Studies may also be feasible for ${\mathchar 28932\relax}^{++}_{{c}{c}}$
baryons.

Experiment | Mode | Yield | Dataset | Ref.
---|---|---|---|---
LHCb | $B_{c}^{+}\rightarrow J/\psi\mu^{+}\nu$ | 19 000 | Run1 | [40]
LHCb | $B_{c}^{+}\rightarrow J/\psi\pi^{+}$ | 25 181 | Run1+Run2 | [41]
LHCb | $B_{c}^{+}\rightarrow J/\psi\pi^{+}\pi^{-}\pi^{+}$ | 9 497 | Run1+Run2 | [41]
LHCb | $B_{c}^{+}\rightarrow J/\psi D_{s}^{+}$ | 1 135 | Run1+Run2 | [41]
LHCb | $B_{c}^{+}\rightarrow B_{s}^{0}\pi^{+}$ | 316 | Run1+Run2 | [41]
CMS | $B_{c}^{+}\rightarrow J/\psi\pi^{+}$ | 7629 | Run2 | [42]
LHCb | ${{\mathchar 28932\relax}^{++}_{{c}{c}}}\rightarrow{{\mathchar 28931\relax}^{+}_{c}}{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}$ | 1598 | Run2 | [43]
LHCb | ${{\mathchar 28932\relax}^{++}_{{c}{c}}}\rightarrow{{\mathchar 28932\relax}^{+}_{c}}{{\pi}^{+}}$ | 616 | Run2 | [43]
Table 6: Yields of doubly-heavy hadrons reconstructed in recent measurements
at the LHC. Run1 corresponds to $\sqrt{s}=7$ TeV and/or 8 TeV, while Run2
corresponds to $\sqrt{s}=13$ TeV.

The predicted fraction, $f_{\rm DPS}\equiv\sigma({{B}_{c}^{+}})_{\rm
DPS}/[\sigma({{B}_{c}^{+}})_{\rm SPS}+\sigma({{B}_{c}^{+}})_{\rm DPS}]$, of
${B}_{c}^{+}$ mesons produced in DPS processes varies as a function of
$p_{\mathrm{T}}$ (Fig. 27), implying the effects would be most pronounced at
low-$p_{\mathrm{T}}$. This would motivate measuring the relative cross
sections as a function of the number of tracks in different $p_{\mathrm{T}}$
regions. The DPS fraction is not found to vary as a function of rapidity.

Figure 27: Fraction of ${B}_{c}^{+}$ decays predicted to be produced by DPS
processes as a function of (left) $p_{\mathrm{T}}$ and (right) rapidity in
simulations samples produced by Pythia.

The contributions from DPS production mechanisms can also be studied in events
with one quarkonium and two singly-heavy hadrons, as discussed in Section 6.
These final states have the advantage that quarkonia can be efficiently
reconstructed using leptonic final states. The singly heavy hadrons could be
reconstructed exclusively, using a range of different final states, or
inclusively, by exploiting the typical topologies of heavy-meson decays.

###  7.1 Exclusive reconstruction of associated singly-heavy hadrons

For singly-heavy hadrons containing $c$-quarks the branching fractions of the
main experimentally-efficient decay channels constitute a reasonable fraction
of the total decay width, as listed in Table 7. Therefore it may be feasible
to exclusively reconstruct a sufficient fraction of the charm hadrons to make
the proposed measurements. Indeed, the cross-section of quarkonia plus one
singly-heavy charm hadron have already been performed [1, 4].

Decay mode | $\mathcal{B}(\%)$
---|---
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ | $3.95\pm 0.03$
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$ | $8.22\pm 0.14$
${{D}^{+}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}$ | $9.36\pm 0.16$
${{D}^{+}_{s}}\\!\rightarrow{{K}^{+}}{{K}^{-}}{{\pi}^{+}}$ | $5.39\pm 0.15$
${{D}^{+}_{s}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$ | $1.08\pm 0.04$
${{\mathchar 28931\relax}^{+}_{c}}\\!\rightarrow{p}{{K}^{-}}{{\pi}^{+}}$ | $6.28\pm 0.32$
Table 7: Branching fractions of experimentally efficient decay channels for
different charm hadron species from Ref [44].

However, for hadrons contain $b$-quarks the typical decay channels have much
smaller branching fractions and may decay further to charm hadrons, therefore
the fraction of events containing a quarkonia in which it is possible to
exclusively reconstruct two $X_{b}$ hadrons is small.

###  7.2 Inclusive reconstruction of associated singly-heavy hadrons

It is possible to inclusively reconstruct singly-heavy hadrons by taking
advantage of the common features of many decay channels, namely a secondary
decay vertex that is significantly displaced from the primary interaction. In
order to measure the kinematic relationships between the singly- and doubly-
heavy hadrons the position of displaced vertices must be measured, along with
a flavour tag labelling the vertex as a $b$\- or $c$-hadron.

Alternatively, the singly-heavy hadrons could be reconstructed as charm and
beauty jets. In this case the reconstruction of the full jet may allow a
better approximation of the parton-level kinematics to be determined rather
than the heavy hadron kinematics. Jet flavour-tagging algorithms use similar
displaced vertex tagging algorithms as already discussed to determine if a jet
is the result of a heavy flavour quark. Recent measurements of heavy dijets at
LHCb claim a reconstruction efficiency for $b\bar{b}$ and $c\bar{c}$ dijets of
about 15% and 1.5%, respectively [45]. Therefore in a significant fraction of
events containing quarkonia it may be feasible to additionally reconstruct
these jets, assuming that the rate of fake jets can be sufficiently
controlled.

##  8 Conclusions

Studies have been performed using the Pythia Monte Carlo event generator to
investigate the contribution from DPS to the production of doubly-heavy
hadrons. New UserHooks have been developed for Pythia to produce events with
one or more heavy quarks more efficiently, enabling the samples to be produced
in a realistic time frame. Comparisons have been made to SPS predictions from
Bcvegpy and GenXicc. The studies show significant and measurable differences
in the ratio of cross sections between doubly- and singly- heavy hadrons as a
function of the collision multiplicity that can discriminate between SPS and
DPS production mechanisms. Further, the DPS mechanism is sensitive to the
assumed scenario for colour reconnections. A range of measurements have been
proposed, including cross-section ratios for ${B}_{c}^{+}$ and ${\mathchar
28932\relax}^{++}_{{c}{c}}$ decays as well as a combination of cross-section
ratios and relative transverse directions in events containing a quarkonium
and two singly-heavy hadrons.

## Acknowledgements

This work is supported by the Monash Warwick Alliance as part of the Monash
Warwick Alliance in Particle Physics (UE, MS, PS, MV), and by the Australian
Research Council via Discovery Projects DP170100708 (PS) and DP210102707 (UE,
TH). MV is supported by the grant ERC-CoG-865469 SPEAR.

## References

  * [1] LHCb collaboration, R. Aaij et al., _Observation of double charm production involving open charm in pp collisions at $\sqrt{s}$ = 7 TeV_, JHEP 06 (2012) 141, arXiv:1205.0975, [Addendum: JHEP 03, 108 (2014)]
  * [2] D0 collaboration, V. M. Abazov et al., _Observation and studies of double $J/\psi$ production at the Tevatron_, Phys. Rev. D90 (2014) 111101, arXiv:1406.2380
  * [3] CMS collaboration, V. Khachatryan et al., _Measurement of prompt $J/\psi$ pair production in pp collisions at $\sqrt{s}$ = 7 Tev_, JHEP 09 (2014) 094, arXiv:1406.0484
  * [4] LHCb collaboration, R. Aaij et al., _Production of associated Y and open charm hadrons in pp collisions at $\sqrt{s}=7$ and 8 TeV via double parton scattering_, JHEP 07 (2016) 052, arXiv:1510.05949
  * [5] ALICE collaboration, B. Abelev et al., _$J/\psi$ production as a function of charged particle multiplicity in $pp$ collisions at $\sqrt{s}=7$ TeV_, Phys. Lett. B712 (2012) 165, arXiv:1202.2816
  * [6] LHCb collaboration, R. Aaij et al., _Study of the doubly charmed tetraquark $T^{+}_{cc}$_, Nature Commun. 13 (2022) 3351, arXiv:2109.01056
  * [7] S. Argyropoulos and T. Sjöstrand, _Effects of color reconnection on $t\bar{t}$ final states at the LHC_, JHEP 11 (2014) 043, arXiv:1407.6653
  * [8] J. R. Christiansen and P. Z. Skands, _String formation beyond leading colour_ , JHEP 08 (2015) 003, arXiv:1505.01681
  * [9] C. Bierlich and J. R. Christiansen, _Effects of color reconnection on hadron flavor observables_ , Phys. Rev. D92 (2015) 094010, arXiv:1507.02091
  * [10] J. Bellm, C. B. Duncan, S. Gieseke, M. Myska, and A. Siódmok, _Spacetime colour reconnection in Herwig 7_ , Eur. Phys. J. C79 (2019) 1003, arXiv:1909.08850
  * [11] T. Sjöstrand, S. Mrenna, and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [12] C. Bierlich et al., _A comprehensive guide to the physics and usage of PYTHIA 8.3_ , arXiv:2203.11601
  * [13] C.-H. Chang, X.-Y. Wang, and X.-G. Wu, _BCVEGPY2.2: A newly upgraded version for hadronic production of the meson $B_{c}$ and its excited states_, Comput. Phys. Commun. 197 (2015) 335, arXiv:1507.05176
  * [14] X.-Y. Wang and X.-G. Wu, _GENXICC2.1: An Improved Version of GENXICC for Hadronic Production of Doubly Heavy Baryons_ , Comput. Phys. Commun. 184 (2013) 1070, arXiv:1210.3458
  * [15] P. Skands, S. Carrazza, and J. Rojo, _Tuning PYTHIA 8.1: the Monash 2013 Tune_ , Eur. Phys. J. C74 (2014) 3024, arXiv:1404.5630
  * [16] T. Sjöstrand and P. Z. Skands, _Transverse-momentum-ordered showers and interleaved multiple interactions_ , Eur. Phys. J. C39 (2005) 129, arXiv:hep-ph/0408302
  * [17] R. Corke and T. Sjöstrand, _Interleaved parton showers and tuning prospects_ , JHEP 03 (2011) 032, arXiv:1011.1759
  * [18] E. Norrbin and T. Sjöstrand, _Production and hadronization of heavy quarks_ , Eur. Phys. J. C17 (2000) 137, arXiv:hep-ph/0005110
  * [19] T. Sjöstrand and P. Z. Skands, _Multiple interactions and the structure of beam remnants_ , JHEP 03 (2004) 053, arXiv:hep-ph/0402078
  * [20] S. Ferreres-Solé and T. Sjöstrand, _The space–time structure of hadronization in the Lund model_ , Eur. Phys. J. C78 (2018) 983, arXiv:1808.04619
  * [21] H. Brooks, C. T. Preuss, and P. Skands, _Sector showers for hadron collisions_ , JHEP 07 (2020) 032, arXiv:2003.00702
  * [22] C.-H. Chang, Y.-Q. Chen, G.-P. Han, and H.-T. Jiang, _On hadronic production of the $B_{c}$ meson_, Phys. Lett. B364 (1995) 78, arXiv:hep-ph/9408242
  * [23] A. V. Berezhnoy, A. K. Likhoded, and M. V. Shevlyagin, _Hadronic production of $B_{c}$ mesons_, Phys. Atom. Nucl. 58 (1995) 672, arXiv:hep-ph/9408284
  * [24] K. Kolodziej, A. Leike, and R. Ruckl, _Production of $B_{c}$ mesons in hadronic collisions_, Phys. Lett. B355 (1995) 337, arXiv:hep-ph/9505298
  * [25] C.-H. Chang, C.-F. Qiao, J.-X. Wang, and X.-G. Wu, _Hadronic production of $B_{c}(B^{*}_{c})$ meson induced by the heavy quarks inside the collision hadrons_, Phys. Rev. D72 (2005) 114009, arXiv:hep-ph/0509040
  * [26] LHCb collaboration, R. Aaij et al., _Measurement of ${B}_{c}^{+}$ production in proton-proton collisions at $\sqrt{s}=$8$\mathrm{\,Te\kern-1.00006ptV}$_, Phys. Rev. Lett. 114 (2015) 132001, arXiv:1411.2943
  * [27] LHCb collaboration, R. Aaij et al., _Measurement of the $B_{c}^{-}$ meson production fraction and asymmetry in 7 and 13 TeV $pp$ collisions_, Phys. Rev. D100 (2019) 112006, arXiv:1910.13404
  * [28] CMS collaboration, _Measurement of production cross sections times branching fraction of $B_{c}^{+}\rightarrow J/\psi\pi^{+}$ and $B^{+}\rightarrow J/\psi K^{+}$ in pp collisions at $\sqrt{s}=7~{}\mathrm{TeV}$ at CMS_, CMS-PAS-BPH-13-002
  * [29] SELEX, M. Mattson et al., _First observation of the doubly charmed baryon $\Xi^{+}_{cc}$_, Phys. Rev. Lett. 89 (2002) 112001, arXiv:hep-ex/0208014
  * [30] LHCb collaboration, R. Aaij et al., _Observation of the doubly charmed baryon ${\mathchar 28932\relax}^{++}_{{c}{c}}$_, Phys. Rev. Lett. 119 (2017) 112001, arXiv:1707.01621
  * [31] T. Sjöstrand and P. Z. Skands, _Baryon number violation and string topologies_ , Nucl. Phys. B659 (2003) 243, arXiv:hep-ph/0212264
  * [32] E598 collaboration, J. J. Aubert et al., _Experimental observation of a heavy particle $J$_, Phys. Rev. Lett. 33 (1974) 1404
  * [33] J.-E. Augustin et al., _Discovery of a narrow resonance in ${e}^{+}{e}^{-}$ annihilation_, Phys. Rev. Lett. 33 (1974) 1406
  * [34] S. W. Herb et al., _Observation of a dimuon resonance at 9.5 GeV in 400-GeV proton-nucleus collisions_ , Phys. Rev. Lett. 39 (1977) 252
  * [35] CDF collaboration, F. Abe et al., _Observation of the $B_{c}$ meson in $p\bar{p}$ collisions at $\sqrt{s}=1.8$ TeV_, Phys. Rev. Lett. 81 (1998) 2432, arXiv:hep-ex/9805034
  * [36] P. L. Cho and A. K. Leibovich, _Color octet quarkonia production_ , Phys. Rev. D53 (1996) 150, arXiv:hep-ph/9505329
  * [37] P. L. Cho and A. K. Leibovich, _Color octet quarkonia production II_ , Phys. Rev. D53 (1996) 6203, arXiv:hep-ph/9511315
  * [38] F. Yuan, C.-F. Qiao, and K.-T. Chao, _$D$ wave heavy quarkonium production in fixed target experiments_, Phys. Rev. D59 (1999) 014009, arXiv:hep-ph/9807329
  * [39] LHCb collaboration, _Evidence for modification of $b$ quark hadronization in high-multiplicity $pp$ collisions at $\sqrt{s}=13$ TeV_, arXiv:2204.13042
  * [40] LHCb collaboration, R. Aaij et al., _Measurement of the ratio of branching fractions $\mathcal{B}(B_{c}^{+}\,\rightarrow\,J/\psi\tau^{+}\nu_{\tau})$/$\mathcal{B}(B_{c}^{+}\,\rightarrow\,J/\psi\mu^{+}\nu_{\mu})$_, Phys. Rev. Lett. 120 (2018) 121801, arXiv:1711.05623
  * [41] LHCb collaboration, R. Aaij et al., _Precision measurement of the $B_{c}^{+}$ meson mass_, JHEP 07 (2020) 123, arXiv:2004.08163
  * [42] CMS collaboration, A. M. Sirunyan et al., _Measurement of B c(2S)+ and B${}_{\mathrm{c}}^{*}$(2S)+ cross section ratios in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, Phys. Rev. D102 (2020) 092007, arXiv:2008.08629
  * [43] LHCb collaboration, R. Aaij et al., _Precision measurement of the $\Xi_{cc}^{++}$ mass_, JHEP 02 (2020) 049, arXiv:1911.08594
  * [44] Particle Data Group, P. A. Zyla et al., _Review of Particle Physics_ , PTEP 2020 (2020) 083C01
  * [45] LHCb collaboration, R. Aaij et al., _Measurement of differential $b\overline{b}$\- and $c\overline{c}$-dijet cross-sections in the forward region of $pp$ collisions at $\sqrt{s}$ = 13 TeV_, JHEP 02 (2021) 023, arXiv:2010.09437