aainstitutetext: Rudolf Peierls Centre for Theoretical Physics, University of
Oxford, OX1 3PU, UKbbinstitutetext: School of Physics and Astronomy, Monash
University, VIC-3800, Australia

# String Junctions Revisited

J. Altmann a,b P. Skands<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recent measurements at the LHC have revealed heavy-flavour baryon fractions
much larger than those observed at LEP, with e.g.,
$\Lambda_{c}^{+}/\mathrm{D}^{0}$ and $\Lambda_{b}^{0}/\mathrm{B}^{0}$ reaching
$\sim 0.5$ at low $p_{\perp}$. One scenario that has been at least partly
successful in predicting observed trends is QCD colour reconnections with
string junctions. In previous work, however, the limit of a low-$p_{\perp}$
heavy quark was not well defined. We reconsider the string equations of motion
for junction systems in this limit, and find that the junction effectively
becomes bound to the heavy quark, a scenario we refer to as a “pearl on a
string”. We extend string-junction fragmentation in Pythia with a dedicated
modelling of this limit for both light- and heavy-quark “pearls”.

##  1 Introduction

The dynamical process by which high-energy quarks and gluons become confined
inside hadrons — _hadronization_ — remains among the most challenging problems
in particle physics. In the context of theoretical models, one typically
constrains a set of non-perturbative parameters in a reference process, like
hadronic $Z$ decays which can be studied cleanly in $e^{+}e^{-}$ collisions,
and then assumes that those same parameters can be reused in different
environments, like $pp$ collisions. This assumption — referred to as _jet
universality_ — is rooted in the factorization of long-distance non-
perturbative physics from short-distance perturbative physics. It underpins,
e.g., the formalism of fragmentation functions, and is also the starting point
for the modelling of hadronization in Monte Carlo (MC) event generators.

It has become clear that there are interesting breakdowns of jet universality
between $e^{+}e^{-}$ and $pp$ collisions and that these breakdowns tend to
become more pronounced with the charged-particle multiplicity of the latter.
Many of the observations, such as enhanced baryon and strange hadron
production CMS:2011fsn ; ATLAS:2011xhu ; CMS:2013zgf ; CDF:2013kip ;
ALICE:2016fzo ; ALICE:2018pal ; ALICE:2023egx , appear reminiscent of
phenomena that are also observed in heavy-ion collisions. Although it is not
yet clear what precise physical conclusions to draw from this, it certainly
motivates a reassessment of the baryon and strange-hadron production
mechanisms in theoretical models, particularly for high-multiplicity $pp$
collisions.

In this work we focus in particular on the Lund string model of hadronization
ANDERSSON198331 ; Andersson:2001yu ; ARTRU1983147 ; Andersson:1998tv , as
implemented in the Pythia 8 event generator Bierlich:2022pfr .

Two modelling aspects that are known to be especially important at high
multiplicities in $pp$ collisions are: multi-parton interactions (MPI) and
colour-space ambiguities / colour reconnections (CR).

A comparatively simple textbook example of the latter, which has been studied
extensively, is hadronic $e^{+}e^{-}\rightarrow WW$ events, in particular in
the context of studying effects of corrections to the so-called Leading-Colour
(LC) string configurations on precision measurements of the $W$ mass
Sjostrand:1993hi ; L3:2003ohc ; L3:2003oci ; OPAL:2003njc ; Siebel:2005uw ;
OPAL:2005cdr ; ALEPH:2006cdc ; DELPHI:2006tie ; ALEPH:2006jcu . However in
$e^{+}e^{-}$ collisions, there is a strong expectation that beyond-LC string
configurations are suppressed, due to several factors including the
$1/N_{C}^{2}$ suppression, the relative boost and space-time separation of the
two $W$ decay systems, and dynamical effects (QCD coherence suppresses “zig-
zagging” colour flows within each $W$ decay system). Nevertheless, the verdict
from a combination of results from all four LEP collaborations was that the
no-CR hypothesis was excluded at 99.5% CL ALEPH:2013dgf .

The picture becomes more complex in $pp$ collision systems (or more generally
hadron-hadron collisions), where one must consider the initial-state coloured
partons, coloured beam remnants, and the contributions of multi-parton
interactions (MPI). As the scattering centres of the MPI sit within $\sim$ a
proton radius of each other, the effect of colour-space ambiguities is
presumably not (significantly) suppressed by space-time separations, nor are
there kinematic or (LC) coherence suppressions. Moreover, there is a
combinatorial enhancement of the ambiguities in events with a large number of
MPI $\sim$ high particle multiplicities, which counteracts the naive
$1/N_{C}^{2}$ suppression.

The very first MC model of MPI for $pp$ collisions already incorporated a
simple CR model Sjostrand:1987su . This was essential to describe the observed
growth of the average charged-particle $p_{\perp}$ with charged-particle
multiplicity, $\left<p_{\perp}\right>(N_{\mathrm{ch}})$, in minimum-bias
collisions. CR also turned out to be crucial for a good description of the
underlying event (UE) at the Tevatron Field:2005sa ; Skands:2010ak . Since
then, many further CR toy models for $pp$ have been proposed Rathsman:1998tp ;
Sandhoff:2005jh ; Buttar:2006zd ; Skands:2007zg ; Gieseke:2012ft ;
Argyropoulos:2014zoa ; Christiansen:2015yqa ; Bellm:2019wrh , often with the
additional motivation to study CR effects on precision determinations of the
top quark mass Skands:2007zg ; Argyropoulos:2014zoa ; Christiansen:2015yqa .
The models that are applied in the $pp$ context are often technically somewhat
simpler than their $ee$ counterparts were, mainly due to the necessity of
addressing more complicated parton topologies. They tend to reconfig. colour
connections based on global potential-energy minimisation arguments. The QCD-
based CR model of Ref. Christiansen:2015yqa  was the first to combine this
with reintroducing the colour-space ambiguities stochastically using SU(3)
colour algebra. It builds randomized SU(3)-weighted colour indices onto the LC
partons from MPI + showers. Thus, full-colour event structures are restored —
if only in an approximate statistical sense. This introduces multiple ways to
achieve colour neutralization, and along with string-length minimization can
allow for alternative string configurations to be selected instead of the LC
ones.

A key feature of the SU(3) structure of QCD is the existence of the
antisymmetric red-green-blue colour-singlet state. As Lund strings form
between colour-connected partons, there naturally must also be a string type
that connects red-green-blue partons in a colour singlet. These strings are
modelled by a Y-shaped structure which we call a “string junction”. The baryon
number of such a junction is conserved such that when junction strings
fragment, baryon formation occurs around the junction itself. Thus the
fragmentation of junction systems provides an additional baryon production
mechanism in $pp$ that is not present in $e^{+}e^{-}$ collision events.
Notably, junction formation scales with the prevalence of CR effects; they are
particularly important for higher-multiplicity events.

In the context of the Lund model, string junctions were first designed to
study baryon-number violating SUSY processes Sjostrand:2002ip , which
correspond to high-$p_{\perp}$ explosions with typically ultra-relativistic
endpoints. The resulting junction-fragmentation formalism was then recycled
for beam remnants, which also involve high-energy endpoints. In neither of
these cases were soft (i.e., low-energy) junction endpoints encountered often
and there was thus no incentive to develop a dedicated treatment of this
limit. However, the same junction modelling was recycled again for QCD CR
Christiansen:2015yqa , but this time due to the explicit string-length
minimization used in the CR procedure, it is fairly frequently the that
topologies are produced in which one or more junction legs are soft. This is
particularly relevant for heavy-flavour baryons at low $p_{\perp}$, which have
received significant experimental interest recently LHCb:2017bdt ;
ALICE:2020wla ; ALICE:2021bli ; ALICE:2021rzj ; LHCb:2021xyh ; ALICE:2022cop ;
ALICE:2022exq ; ALICE:2023jgm ; ALICE:2023xiu ; ALICE:2023wbx ; ALICE:2023sgl
; LHCb:2023wbo . The aim of this work is therefore to consider this limit more
carefully than was done in the past, enabling us to make firmer predictions.

We proceed systematically, and consider low-energy excitations around the
junction in general, including also how soft gluon kinks affect the effective
juntion motion and the resulting junction baryon spectra. This paper is
organised as follows. In sec. 2 we summarize briefly the Lund String Model
with particular focus on string junctions. In sec. 3 we elaborate on the
modelling of junction motion, exploring both the existing modelling and
previously unconsidered so-called soft-leg cases. We also provide a model of
gluon kinks passing through junctions, however this is limited to a
theoretical description and is not implemented in Pythia. In sec. 4 we
describe the updated implementation in Pythia, addressing the treatment of
soft gluon kinks on junction motion. Finally, sec. 5 outlines some comparison
of the updated implementation to both the old junction modelling and to data.

##  2 Background

###  2.1 Lund String Model

The Lund String Model is a semi-classical phenomenological hadronization
model, based on confinement in Quantum Chromodynamics (QCD). The model
collapses the colour confinement field into an infinitely narrow flux tube
spanned between coloured particles, modelled as a 1+1 dimensional relativistic
worldsheet which we call a string, and is modelled using relativistic string
dynamics ANDERSSON198331 ; Andersson:2001yu ; ARTRU1983147 . Modelling the
confinement field as a string is largely motivated by the Cornell potential
Eichten:1978tg ; Bali:1992ab , a lattice QCD result that demonstrates a linear
potential between static colour charges. Such a linear potential corresponds
to a constant tension, and hence a string characterised by the constant string
tension $\kappa\sim 1$ GeV/fm.

These confinement field form between colour-connected partons, such that the
overall string configuration is in a singlet state according to the SU(3)
colour structure of QCD. This means that quarks/antiquarks (i.e.
triplet/antitriplet colour charges) are colour connected to a single other
parton, and gluons (i.e. octet colour charges) will be colour connected to two
other partons via two string pieces. The simplest string form to consider is
the so-called dipole string, which makes use of the colour-anticolour colour
neutral state and has a quark and antiqaurk endpoint, with arbitrarily many
intermediate gluons that form transverse “kinks” on the string. An example of
such a dipole string with two gluon kinks is shown in fig. 1.

Figure 1: Sketch of a string configuration with a single gluon-kink on a
dipole string. Here the strings are represented by the dashed coloured lines
and are shown at several different times. Here one can see that a gluon-kink
results in two separate string segments, one connecting $g_{R\bar{B}}$ to
$q_{B}$ and another connecting $g_{R\bar{B}}$ to $\bar{q}_{\bar{R}}$.

####  2.1.1 Beyond Leading Colour

In event generators like Pythia, the mapping of colour flow is simplified by
considering the Leading-Colour (LC) limit, in which QCD is modelled as an
$SU(N_{C})$ gauge theory with the number of colours taken to infinity (i.e.
$N_{C}\rightarrow\infty$). In this limit, gluon (adjoint) colour states can be
reduced to simple direct products of fundamental and antifundamental ones,
since the singlet can be neglected in the group relation

$N_{C}\otimes\bar{N}_{C}~{}=~{}\overbrace{(N_{C}^{2}-1)}^{\mathrm{adjoint}}~{}\oplus~{}1~{}.$
(1)

This vastly simplifies the problem of colour flow in parton cascades and also
eliminates colour interference effects beyond the dipole level; each colour in
an event is unique and is matched to a single other anticolour. At the
perturbative level, these colour connections are used to set up the system of
radiating LC “dipoles”, while at the non-perturbative level the LC colour
connections dictate between which partons confining potentials should arise:
in the Lund string model, each LC dipole is dual to a non-perturbative string
piece Gustafson:1986db . Thus, the LC limit provides for unambiguous string
topologies, with each string piece having a colour and an anticolour charged
endpoint, i.e. dipole strings.

When generalising from a single parton system to multi-parton interactions,
the LC limit does not address if or how different MPI systems should be
correlated in colour space. In the limit of high $Q^{2}$, one may assume that
each MPI is independent of any others in colour space. In an LC picture, this
would imply that there would be no string connections directly between
outgoing partons from different (high-$Q^{2}$) multi-parton interactions
(MPIs). (At lower $Q$ of order the inverse proton size, coherence and/or
saturation effects presumably modify this picture, but not much is known about
the details.)

What about effects beyond LC? In hadronic $Z$ decays, effects beyond LC are
expected to be very small, as they are suppressed both by $1/N_{C}^{2}\sim
10\%$ and further by a combination of kinematics and coherence (no overlapping
jets and dominance of angular ordered colour structures within each jet).
However, in systems with several independent colour sources, such as in
hadronic $e^{+}e^{-}\to WW$ events (in the limit of vanishing $W$ lifetime111
Presumably a reasonable starting approximation since $\Gamma_{W}\sim
2\,\mathrm{GeV}\,\gg\,\Lambda_{\mathrm{QCD}}$.) the kinematic suppression can
be lifted in parts of phase space, and in dense string environments such as
(high-multiplicity) $pp$ collisions, combinatorics can further counteract the
naive $1/N_{C}^{2}$ suppression. That is, while effects beyond LC between two
dipoles, say, are suppressed by $1/N_{C}^{2}$, the probability that there are
no beyond-LC effects in a system of $n$ dipoles will scale like
$(1-1/N_{C}^{2})^{n-1}$, which becomes asymptotically small for large $n$.

Focusing now on QCD, with $N_{C}=3$, the finite number of colours allows for a
potentially large number of _possible_ string configurations to produce
overall singlet states, resulting in ambiguity in where the confining fields
form. The QCD-based colour reconnection (CR) model of Ref.
Christiansen:2015yqa  stochastically reintroduces such colour-space
ambiguities by assigning colour/anticolour indices from 0 to 8 to partons in
order to reproduce probabilities defined by the $\mathrm{SU}(3)$ algebra,

$3\otimes\bar{3}=8\oplus 1~{},$ (2a) $3\otimes 3=6\oplus\bar{3}~{},$ (2b)

where the octet and sextet states are interpreted as non-confining, while the
singlet and antitriplet are interpreted as confining and partially confining
respectively.

To determine which among the many resulting possible string configurations is
realised in terms of being the one selected for hadronisation, the
configuration that gives the overall “shortest string lengths” is chosen. Here
the “length” of the string actually refers to a momentum-space Lorentz-
invariant measure of the integrated energy density per unit length of the
string, which we call the $\lambda$-measure. In the context of the QCD CR
model, the default form for the $\lambda-$measure for a $q\bar{q}$ dipole
string was previously

$\lambda^{q\bar{q}}=\ln{\left(1+\frac{\sqrt{2}E_{1}}{m_{0}}\right)}+\ln{\left(1+\frac{\sqrt{2}E_{2}}{m_{0}}\right)}~{},$
(3)

where $E_{i}$ is the energy of the parton in the dipole rest frame, and
$m_{0}$ is a regularisation parameter of order $\Lambda_{\mathrm{QCD}}$. For
massive endpoint partons, however, this form overestimates the physical string
length, since the endpoint masses are allowed to contribute to the energies in
the numerators.

To ensure a more sensible treatment of QCD CR involving strings with heavy-
quark endpoints, we introduce the following generalisation (made default from
Pythia 8.311 onwards),

$\lambda^{q\bar{q}}=\text{max}\left[\ln{\left(\frac{E_{1}+|\vec{p}|}{m_{1}+m_{0}}\right)},0\right]+\text{max}\left[\ln{\left(\frac{E_{2}+|\vec{p}|}{m_{2}+m_{0}}\right)},0\right]~{},$
(4)

with $m_{i}$ the mass and $(E_{i},\pm\vec{p})$ the energies and momenta of the
partons in the rest frame of the pair. This measure ensures the limit of
$|\vec{p}|\rightarrow 0$ produces $\lambda\rightarrow 0$. Here $m_{0}$ is used
to protect against massless partons, and the max functions impose that no
endpoint can be associated with a negative contribution to the effective
$\lambda$ measure.

We note that an alternative $\lambda$-measure for string pieces involving
heavy-quark endpoints was already proposed in Bierlich:2023okq , based on the
rapidity of a heavy-quark endpoint as a measure of its associated string
length. But that measure only applied to heavy quarks and could not be used
for light or massless partons. The measure defined in eq. (4) is similar in
spirit to the one proposed in Bierlich:2023okq  but with modifications so that
the same form of the $\lambda$-measure can be used for all partons regardless
of mass, allowing consistency across all types of string lengths being
compared in the CR treatment. This same string-length measure is then also
used to calculate lengths for $qg$, $gg$ and $g\bar{q}$ string connections as
well.

By allowing colour reconnections, string connections between different MPI can
be formed. Moreover, in addition to strings spanned between triplet and
antitriplet endpoints, two other string topologies can arise: gluon loops
Andersson:1984af  and junctions Sjostrand:2002ip , illustrated in fig. 2.

Figure 2: Types of string topologies, with “textbook examples” of each kind.

The latter will be the focus of this paper and revisited in sec. 2.2,
including details on the $\lambda-$measure for junction topologies and the
fragmentation of such strings.

####  2.1.2 Fragmentation

In high-energy collisions such as those at the Large Hadron Collider (LHC),
colour-connected partons move apart at high energies, resulting in high-
invariant-mass string systems that stretch the confining potentials to the
point of breaking. At the site of these string breaks, a quark-antiquark (or
diquark-antidiquark) pair is created. This string-breaking process is called
fragmentation, and occurs along the string until there is no longer sufficient
energy to keep fragmenting the string, resulting in final-state “primary”
hadrons (some of which may be unstable and undergo further decays, creating
“secondary” hadrons). A basic schematic of this process is shown in Fig 3.

Figure 3: A basic schematic showing the hadronization process at different
times for a simple dipole string. i) The initial colour-anticolour pair begin
to move apart at high energies. ii) The colour flux tube forms between the
colour-anticolour charge pair. iii)-iv) New $q\bar{q}$ pairs are formed from
string breaks. v) After hadronization is complete, we are left with colourless
hadrons which form jets.

The modelling of string fragmentation in Pythia (the Lund model) relies on two
main components; the quantum tunnelling process for spontaneous pair creation,
which dictates the flavour and transverse components of the fragmentation
process, and the fragmentation function that governs the longitudinal
component. To model the spontaneous pair creation that occurs in a string
break, a QCD analogy of the QED Schwinger mechanism PhysRev.82.664  is used.
The Schwinger mechanism was originally derived in the context of spontaneous
electron-positron pair creation in the presence of a strong electric field. In
the Lund model, an analogous QCD formulation is used for $q\bar{q}$ creation
from the confinement field. The leading term of the Schwinger mechanism
implies a Gaussian form with respect to the transverse mass $m_{\perp}$,

$\exp\left(\frac{-\pi m_{\perp}^{2}}{\kappa}\right)=\exp\left(\frac{-\pi
p_{\perp}^{2}}{\kappa}\right)\exp\left(\frac{-\pi m^{2}}{\kappa}\right),$ (5)

with parton mass $m$ and transverse momentum $p_{\perp}$ relative to the local
string axis. This provides a suppression of both heavy-particle production and
constrains the $p_{\perp}$ distribution. Alternative string-breaking models
have also been proposed (e.g., the thermal model of Ref. Fischer:2016zzs ),
however so far the standard has remained Schwinger-type string breaks. This is
what will be used throughout the remainder of this paper. An important
consequence of the Gaussian suppression and the string tension $\kappa\sim 1$
GeV/fm is that only light-flavoured quarks (i.e. up, down and strange quarks)
can be created via fragmentation. Thus all charm and beauty quarks must come
from perturbative processes (which can include MPI), a distinction of
particular importance as it allows heavy-flavour hadrons to serve as
interesting probes of fragmentation modelling. We will return to this point in
sections 2.2 and 3.

The longitudinal component of fragmentation is governed by what is known as
the Lund symmetric fragmentation function (or also the left-right symmetric
fragmentation function). In the Lund string model, string breaks are treated
as causally disconnected, and hence the time ordering of string breaks holds
no physical significance. This allows us to perform string breaks in any order
we choose. The simplest is to fragment hadrons off from either string
endpoint, which makes it particularly easy to impose on-shell hadron mass
constraints on each produced hadron. Formulated with a simple dipole string in
mind, the symmetric fragmentation function takes the form

$f(z)=N\frac{1}{z}(1-z)^{a}\exp{\left(\frac{-bm_{\perp}^{2}}{z}\right)},$ (6)

which is a probability distribution for the hadron to take fraction $z$ of the
string momentum, with tuneable free parameters $a$ and $b$, and $N$ a
normalisation constant. The string end from which to fragment is normally
chosen randomly for each break.

###  2.2 Junction String Topologies

String junctions arise naturally from the SU(3) colour structure of QCD and
represent the confinement field spanned between three colour-triplet partons
(or three colour-antitriplet ones) in an overall colour-neutral state. These
states only exist for finite $N_{C}$ and hence intrinsically go beyond the LC
picture. In the context of the Lund model, with $N_{C}=3$, they are
represented by Y-shaped string topologies Sjostrand:2002ip  as illustrated in
the right-hand side of fig. 4.

Figure 4: The left image shows an LC string configuration consisting of three
dipole strings. The coloured lines here represent the strings, with arrows
indicating the direction of the colour flow (conventionally flowing from
colour to anticolour). The right side image shows a possible alternative
string configuration given a junction-type colour reconnection, resulting in
the formation of a junction and an antijunction string system.

Here the “junction” itself refers to the topological feature depicted by the
central vertex point of the Y-shape, with the “junction legs” being the three
string pieces extending out from the junction.

In order for string junctions to be formed via colour reconnections, one needs
a generalised string-length measure which allows to compare the effective
string lengths between dipole- and junction-type configurations, such as those
on the left and right of fig. 4 respectively. With the measure for dipole-type
string pieces defined by eq. (4) above, we need an equivalent measure for
parton-to-junction string pieces. This is constructed by reinterpreting each
junction leg as half of a dipole string, with the midpoint of the dipole
sitting at the junction itself (in a frame in which the junction is at rest).
In a straightforward generalisation of eq. (4), we define the
$\lambda$-measure associated with each of the parton–junction string segments
in a three-parton junction configuration, as

$\lambda^{qqq}_{i}=\text{max}\left[\ln{\left(\frac{E_{i}+|\vec{p}_{i}|}{j_{0}(m_{i}+m_{0})}\right)},0\right]~{},$
(7)

Importantly, the energies $E_{i}$ and momenta $\vec{p}_{i}$ of each of the
three partons are defined in the rest frame of the junction, which generalises
the notion of the dipole rest frame that was used in eq. (4). What this
junction rest frame (JRF) looks like and how to find it will be explored in
the detail in sec. 3.

The $\lambda$-measure defined in eq. (7) also introduces a free parameter,
$j_{0}$, which allows the user to modify how easily junction systems form.
Larger values of $j_{0}$ make junction formation more likely (by decreasing
the effective $\lambda$-measures for junction topologies relative to dipole
ones), and vice versa. In the code, $j_{0}$ is set by the tuneable parameter
ColourReconnection:junctionCorrection. As an option mainly intended for
comparisons, we also retain the possibility to choose to use the old form of
the string-length measure instead, as per eq. (3), where this parameter
$j_{0}$ multiplies the $m_{0}$ term in the denominator.

###  2.3 Junction Fragmentation

The string-length measures defined above govern whether and how easily
junction topologies are formed in CR. The next question is how does one
fragment a junction string system? The standard approach to junction
fragmentation Sjostrand:2002ip  is to use the concept of reinterpreting
junction legs as equivalent to half of a dipole string, analogously to how the
string-length measure was defined above. This allows us to simply recycle the
dipole string-fragmentation procedure from sec. 2.1.2. Construction of the
full dipole string is done by boosting the junction system to the junction
rest frame, and creating what we call a fictitious leg. This fictitious leg is
a mirror image of a junction leg that extends on the opposite side of the
junction, a depiction of which can be seen in fig. 5. This fictitious leg acts
as the combined effective pull of the other two junction legs, providing a
reservoir of oppositely oriented momentum.

Of course we generally encounter more complicated junction systems than those
shown in fig. 5, which can include one or more gluon kinks on each of the
junction legs. In such cases, the rest frame of the junction will change over
time as different gluon-kinks dominate the junction motion. Rather than
fragmenting the junction legs with a dynamically changing junction velocity,
an average JRF is determined and used in order to construct the fictitious
legs. However, the original procedure for determining the average JRF proposed
in Ref. Sjostrand:2002ip  often failed to converge, especially when (slow)
heavy quarks were involved. In this work, we develop a new and more stable
procedure for determining the average JRF, described in sec. 3. However note
that this average JRF is only used for fragmentation and not when deciding
where colour reconnections are formed. (For CR string-length calculation
purposes, only the “first frame”, defined by the partons nearest to the
junction, is considered.)

Figure 5: Fragmentation of a junction topology using the standard junction
fragmentation procedure in the Mercedes frame. The first frame shows the
fragmentation of the two legs with the lowest momenta in the JRF (which here
is the Mercedes frame). For these junction legs, fictitious legs are
constructed on the opposite side of the junction. These are shown by the
dashed lines with endpoints $\bar{q}_{01}$ and $\bar{q}_{02}$. These two legs
have been fragmented from their real endpoints inwards towards the junction,
after which we have the formation of the junction diquark, $q_{3}q_{5}$, from
the string breaks next to the junction. The right image shows fragmentation of
the final junction leg, treated as a dipole with endpoints $q_{2}q_{5}$ and
$q_{03}$. Here the strings breaks next to the junction, $q_{2}q_{5}q_{9}$
combine to form the so-called junction baryon.

Once the average JRF has been found and the fictitious endpoints constructed,
the two softest junction legs are fragmented in towards the junction, cf. the
left side of fig. 5. In Pythia 8.310 and prior, the two softest junction legs
are determined by the legs with the lowest JRF energies. In this work and from
Pythia 8.311 onwards, the lowest absolute momentum is used instead. This
change was made as the energy of a junction leg does not reflect how soft the
leg is for heavy-quark endpoints.

In the example depicted in fig. 5, the two softest junction legs are defined
by endpoints $q_{01}$ and $q_{02}$. Once these junction legs are identified
and their fictitious endpoints constructed, each of these junction legs is
fragmented from the real endpoint in towards the junction, using standard
string-fragmentation methods (except for not alternating between the two
ends), with fragmentation stopping once the junction is reached. To measure
when to stop fragmentation of a given leg, the energy of the leg in the JRF is
compared to the summed energies of the produced hadrons. There are also
parameters StringFragmentation:eBothLeftJunction and
StringFragmentation:eMaxLeftJunction which regulate the remaining energy of
the two softest junction legs in the JRF so that it is not too large.

Additional to this energy measure, as of Pythia 8.311 we have also included a
further constraint which ensures that the mass of the so-called “junction
diquark” remains positive. This junction diquark is comprised of the quarks
from the string breaks closest to the junction from the fragmentation of the
two softest junction legs. In fig. 5, this junction diquark is defined by
$q_{2}q_{5}$. For the momentum of this junction diquark, energy-momentum
conservation yields

$p_{\text{diq}}=p_{\text{min}}+p_{\text{mid}}-p_{\text{had}},$ (8)

where $p_{\text{min}}$ and $p_{\text{mid}}$ are the momenta of the softest and
next-to-softest junction legs respectively, and $p_{\text{had}}$ is the
momentum sum of the produced hadrons from the fragmentation of these two legs.
Once the junction diquark is formed, the last junction leg can be hadronized
as a standard dipole string, which in the example in fig. 5 would have
endpoints $q_{3}q_{5}$ and $\bar{q}\bar{q}_{03}$. Note that the choice of
fragmenting the softest two junction legs first is somewhat arbitrary, however
allows easier handling of hadron mass constraints as it avoids the dipole
constructed from the junction diquark and final junction endpoint from being
too small. The impact of making alternative ordering choices was studied in
Ref. Sjostrand:2002ip  and found to be small.

Notably one can see in fig. 5 that the string breaks surrounding the junction
necessarily form a baryon, which here we call the “junction baryon”. This
provides an additional baryon production mechanism unique to beyond-LC
modelling of colour reconnections plus string fragmentation (with standard
baryon production resulting from diquark-antidiquark string breaks or beam
remnants). The significant effect of junction baryon production has already
been seen in increased baryon-to-meson ratios in $pp$ collisions. An early
example of this was the $\Lambda/K$ ratio measured by CMS CMS:2011jlm , which
— as shown in Ref. Christiansen:2015yqa  — default Pythia (which relies on
LEP-LHC universality Skands:2014pea ) underpredicts by over 20% while the CR
model with junctions fits the data well.

Though junction modelling is relevant for all baryon production, it can be
particularly relevant for heavy-flavour baryon production. Heavy-flavour
(charm and beauty) quarks act as interesting probes into hadronization as
heavy flavours cannot themselves be created via string breaks. Thus in order
to form a heavy-flavour junction baryon, the junction leg containing the
heavy-quark endpoint must be sufficiently soft that it does not fragment but
remains part of the junction baryon. This leaves the heavy-flavour baryon
sensitive to the modelling of the junction motion and junction fragmentation.

Large effects of junctions on baryon production in the heavy-flavour sector
have already been seen, particularly at low $p_{\perp}$. This effect can be
seen in the increased $\Lambda_{c}/D^{0}$ ratio at low $p_{\perp}$
ALICE:2021rzj , where the default Pythia tune (Monash 2013 Skands:2014pea )
predicts an approximately flat distribution. Given that the probability of
diquark-antidiquark pair production via string breaks has some fixed
probability (in the absence of collective effects such as ropes
Bierlich:2016faw ; Bierlich:2017sxk  or close-packing Fischer:2016zzs ), the
universality-limit flat distribution is unsurprising. The rise at low
$p_{\perp}$ exhibited by the QCD CR models can be explained by junction
baryons as junctions predominantly sit at low $p_{\perp}$. This is a natural
consequence of the string-length minimisation in CR, which results in junction
reconnections largely occurring between jets the tips of which fragment into
hard mesons (and baryons with the normal — universal — ratio) while the
junction baryon is the most subleading and hence typically softest hadron
produced (in the JRF), hence it generally also sits at low $p_{\perp}$ (in the
LAB).

Another interesting consequence of including junctions in a string model is
the decrease of the $\Lambda_{b}$-baryon asymmetry at low $p_{\perp}$ as seen
by LHCb LHCb:2021xyh . For $pp$ collisions, models without junction CR predict
a large asymmetry. This arises from $\Lambda_{b}$ production via combination
of a $b$ quark with the proton beam remnant, whereas there is no equivalent
mechanism for $\bar{\Lambda}_{b}$ due to the absence of an antiproton beam
remnant for a $\bar{b}$ quark to combine with. However in CR models with
junctions, junctions are always created in pairs (total baryon number is
conserved), cf., e.g., fig. 4. The contributions from junction topologies to
$\Lambda_{b}$ and $\bar{\Lambda}_{b}$ production are therefore approximately
equal, and as previously discussed particularly occurs at low $p_{\perp}$
which in turn dilutes the asymmetry in that part of phase space.

The original formulation of junction fragmentation, however, focused on high-
energy junction legs Sjostrand:2002ip . Special consideration of the limit of
soft endpoints was not built into the model, nor were mass effects arising
from gluon kinks and endpoint masses under good control. Consequently, this
initial construction, which remained the implementation for Pythia 8.310 and
earlier, was not stable against such cases, with the procedure in Pythia 8.310
unable to construct an average JRF for around 10% of minimum-bias events at
LHC energies. The interesting predictions seen for low-$p_{\perp}$ heavy-
baryon production with the QCD CR model including junctions motivates careful
examination of such soft endpoint junctions. Also given that CR aims to
minimise string lengths, one would naturally expect many junction systems to
contain soft endpoints and soft gluon kinks near the junction. Thus in order
for us to make firmer predictions, particularly but not only for
low-$p_{\perp}$ heavy baryons, we here aim to construct a more physically
robust model. The first step is to reexamine the determination of the average
JRF.

##  3 Junction Motion

The most intuitive way to construct the JRF is to consider balancing the force
exerted by the three junction strings on each other. Given the constant force
from strings and a simple three-parton junction configuration, the natural
frame to consider is where the opening angle between the 3-momenta of each
pair of partons is 120° Sjostrand:2002ip . In the following we will refer to
this as the “Mercedes frame”, an example of which can be seen in fig. 5. For a
set of massless four-vectors, a boost to the Mercedes frame always exists (we
will elaborate on this point in sec. 3.1), however the same is not true for
massive four-vectors. In the limit the Mercedes frame no longer exists, by
considering the classical action of the string, the junction effectively
becomes bound to the massive parton and follows its motion. This scenario,
which we call “pearl-on-a-string”, will be described in detail in sec. 3.2.2.

In the following, we begin by describing the Mercedes frame and then proceed
to consider soft-leg cases, examining both oscillatory motion of an endpoint
around a junction and the pearl-on-a-string scenario. In these discussions, we
will begin by remaining in the context of a simple three-parton junction
configuration without the inclusion of gluon kinks. Then, we introduce a
single gluon kink on a junction leg in sec. 3.3 and provide a theoretical
description of the effect on junction motion. The practical handling of
junction configurations generalised to multiple gluon-kinks is explored in
sec. 4.

###  3.1 Mercedes frame

The boost to the Mercedes frame can be determined using the method presented
in Sjostrand:2002ip , which is summarised for context below (note this
procedure is unchanged in the updated implementation). Let us label the parton
four-vectors $p_{i}$ in the initial (arbitrary) frame of reference, with
$p^{\prime}_{i}$ being the four-vectors in the Mercedes frame. Given the
Lorentz invariant four-products, $a_{ij}=p_{i}p_{j}$, and the fixed angle
between the three-momenta in the Mercedes frame,
$\theta_{ij}=120^{\circ}=2\pi/3$, the Lorentz invariant four-products can be
defined in terms of the Mercedes-frame energy and momentum,

$a_{ij}=p_{i}p_{j}=p^{\prime}_{i}p^{\prime}_{j}=E^{\prime}_{i}E^{\prime}_{j}-|\vec{p\,}^{\prime}_{i}||\vec{p\,}^{\prime}_{j}|\cos{\frac{2\pi}{3}}=E^{\prime}_{i}E^{\prime}_{j}+\frac{1}{2}|\vec{p\,}^{\prime}_{i}||\vec{p\,}^{\prime}_{j}|.$
(9)

By rewriting the energies in terms of the mass and momenta, one can see that
the only degrees of freedom in eq. (9) are $|\vec{p\,}^{\prime}_{i}|$ and
$|\vec{p\,}^{\prime}_{j}|$. Thus by rearranging eq. (9) we can introduce a
function, $f_{ij}$, such that solutions for $|\vec{p\,}^{\prime}_{i}|$ and
$|\vec{p\,}^{\prime}_{j}|$ are found when $f_{ij}=0$. Note that the first two
arguments of this function are monotonically increasing.

$f_{ij}=f(|\vec{p\,}^{\prime}_{i}|,|\vec{p\,}^{\prime}_{j}|;m_{i},m_{j},a_{ij})=\sqrt{|\vec{p\,}^{\prime}_{i}|^{2}+m_{i}^{2}}\sqrt{|\vec{p\,}^{\prime}_{j}|^{2}+m_{j}^{2}}+\frac{1}{2}|\vec{p\,}^{\prime}_{i}||\vec{p\,}^{\prime}_{j}|-a_{ij}$
(10)

By requiring $f_{12}=0$ and $f_{13}=0$ whilst allowing
$|\vec{p\,}^{\prime}_{1}|$ to vary freely in a kinematically allowed region,
we can uniquely solve for the other momenta as a function of
$|\vec{p\,}^{\prime}_{1}|$,

$|\vec{p\,}^{\prime}_{j}|(\vec{p\,}^{\prime}_{1})=\frac{2E^{\prime}_{1}\sqrt{4a_{1j}^{2}-m_{j}^{2}(4E^{{}^{\prime}2}_{1}-|\vec{p\,}^{\prime}_{1}|^{2})}-2|\vec{p\,}^{\prime}_{1}|a_{ij}}{4E^{{}^{\prime}2}_{1}-|\vec{p\,}^{\prime}_{1}|^{2}}.$
(11)

Using eq. (11), we can rewrite
$f_{23}(|\vec{p\,}^{\prime}_{2}|,|\vec{p\,}^{\prime}_{3}|)$ as a function of
$|\vec{p\,}^{\prime}_{1}|$. As eq. (11) is a decreasing function with
$|\vec{p\,}^{\prime}_{1}|$, the function $f_{23}$ is also a monotonically
decreasing function with $|\vec{p\,}^{\prime}_{1}|$. Thus we can use
$f_{23}(|\vec{p\,}^{\prime}_{1}|)=0$ to solve for $|\vec{p\,}^{\prime}_{1}|$
and therefore solve for all $|\vec{p\,}^{\prime}_{i}|$. In the case where all
three junction legs are massless, there is always a physical solution which is
given by

$E^{\prime}_{i}=|\vec{p\,}^{\prime}_{i}|=\sqrt{2a_{ij}a_{ik}/3a_{jk}}~{}.$
(12)

For general masses, however, there is no simple analytical solution and
instead iterative root-finding procedures are used. But for some three-parton
configurations involving at least one massive parton, a physical solution to
the above procedure may not exist at all, implying that a boost to the
Mercedes frame does not exist. That is, for certain parts of massive three-
body phase space, there is no Lorentz frame in which the opening angles
between all three partons are 120 degrees. These cases constitute one class of
“soft-leg" cases, and will be treated separately (see below in sec 3.2.2).
Note that in Pythia 8.310 and prior, there is no special consideration for
such scenarios, and the root-finding procedure would generally terminate in an
error.

When the Mercedes frame exists and having determined the energies and momenta
of each parton in that frame according to the above procedure, one can easily
construct a boost to this frame. We use the centre-of-mass frame as a stepping
stone as it is a well-defined frame and the momenta in this frame lie in a
plane, which necessarily will be the same plane the Mercedes frame sits on.
The boost to the centre-of-mass frame is given by
$\vec{\beta\,}^{CM}=-\sum\vec{p}_{i}/\sum E_{i}$, with centre-of-mass momenta
$p_{i}^{CM}$. Thus for boost $\vec{\beta\,}^{\prime}$ from the centre-of-mass
frame to the Mercedes frame, $\vec{\beta\,}^{\prime}$ must obey
$\gamma^{\prime}E_{i}^{CM}+\gamma^{\prime}\vec{\beta\,}^{\prime}\vec{p\,}_{i}^{CM}=E^{\prime}_{i}$.
Dividing this equation by $E_{i}^{CM}$, and subtracting the equation for $i=j$
from $i=1$, we get

$\gamma^{\prime}\vec{\beta\,}^{\prime}\biggl{(}\frac{\vec{p\,}_{1}^{CM}}{E_{1}^{CM}}-\frac{\vec{p\,}_{j}^{CM}}{E_{j}^{CM}}\biggr{)}=\frac{E^{\prime}_{1}}{E_{1}^{CM}}-\frac{E^{\prime}_{j}}{E_{j}^{CM}}.$
(13)

From here we can parameterize the boost to the Mercedes frame as a linear sum
of the vector differences in eq. (13) for $j=2$ and $j=3$. Then combining
$\beta^{CM}$ and $\beta^{\prime}$ gives the overall boost to the Mercedes
frame from the original frame.

###  3.2 Soft-leg Junction Motion

We define a junction leg as “soft” if it does not have sufficient energy for a
string break to occur between the endpoint and the junction. In the following,
we will consider junction motion in the absence of such a break, and will
assume there is only a single soft leg. We consider both massless and massive
limits for the soft leg, while we shall assume the other two legs to have
massless endpoints and be much more energetic in comparison.

####  3.2.1 Massless Soft-leg Case

Consider the case of a single soft massless endpoint in the Mercedes frame. In
the absence of string breaks, it will travel outwards from the junction until
it has lost all of its momentum to the junction string. After that, it will
change its direction of motion and begin moving back towards the junction.
Eventually, it will “hit” the junction, at which time the junction will no
longer remain at rest and instead begins moving in the direction of the soft
endpoint. This process is repeated when the parton turns around again. Overall
this behaviour results in some oscillation of the soft leg around the
junction, with the junction itself moving in a sort of start-stop motion, at
rest when the soft parton is on one side of it and moving when it is on the
other.

To study this oscillatory motion in detail, it is instructive to examine the
behaviour in the so-called Ariadne frame with respect to the soft junction
leg. As depicted in the first frame of fig. 6, the Ariadne frame is defined
such that the 3-momenta of the two more energetic legs are back-to-back (here
we allign this with the $z$-axis) and the soft leg is at 90° to the other two
legs (along the labelled $x$-axis in this case). This is a convenient frame to
map the junction motion as the $x$-direction contribution is solely due to the
soft endpoint.

Figure 6: Given a simple three-parton junction configuration in the Ariadne
frame with respect to $q_{1}$, here we show the oscillatory behaviour of
parton $q_{1}$ around the junction at different times. This sketch is drawn
assuming no string breaks occur and all massless partons. Here the initial
momentum of $q_{1}$ is $p_{0}$, and the distance it traverses till it loses
its momentum in the first iteration is given by $t=p_{0}/\kappa$. The initial
momenta of $q_{2}$ and $q_{3}$ are assumed to be much greater than $p_{0}$.
Note that the arrows here are only intended to illustrate the momentum
directions, not their magnitudes.

Given that the initial momenta of each parton are massless for now, the
Mercedes frame exists and defines the JRF at early times. From a Mercedes
frame supposing massless partons, a boost of $\beta=1/2$ is required to get to
the Ariadne frame, implying a junction velocity of $1/2$ in this frame. Using
the initial junction velocity of 1/2 and parton velocity of 1, the oscillatory
motion can be fairly simply mapped out as in fig. 6.

In fig. 6, the initial 3-momentum of the soft endpoint $q_{1}$ is aligned with
the positive $x$ direction, $\vec{p}_{0}=p_{0}(1,0,0)$. The equation of motion
for a string endpoint is $\mathrm{d}p/\mathrm{d}t=-\kappa$ Andersson:1998tv
(oriented along the string). Thus, the time it takes for $q_{1}$ to lose all
of its momentum, and hence the distance it will traverse before it must turn
around, is linearly proportional to its momentum. We label this time as
$t=p_{0}/\kappa$.

The second image in fig. 6 shows the configuration at time $t$, at which
$q_{1}$ has lost its momentum and will begin to move in the negative
$x$-direction. This change has not propagated to the junction yet, however,
which will continue to move in the positive $x$-direction with velocity
$\beta=1/2$. With these new velocities, one can calculate that after a further
time $t/3$ the junction will catch up to parton $q_{1}$ as seen in the third
image of fig. 6. By this time $q_{1}$ will have gained a momentum
$-\vec{p}_{0}/3$ from the string.

As $q_{1}$ passes through the junction, the net force exerted on the junction
from junction strings will change, and the junction will now begin to move in
the negative $x$-direction with velocity $\beta=-1/2$. This forms a zig-zag
shape on the string around the junction as seen in fourth image of fig 6. One
can see that the shape of strings close to the junction after the zig-zag has
formed is the same shape as at initial junction string configuration, however
scaled down by a factor of $1/3$.

From here, the same calculations can be used to map out the next oscillation
as shown in the bottom three frames of fig 6, whereby given momentum
$-\vec{p}_{0}/3$ at the start of the oscillation, $q_{1}$ will now lose its
momentum after a further time $t/9$ and will next hit the junction with
momentum $\vec{p}_{0}/9$. This oscillatory pattern then continues, scaled down
by a factor of $1/3$ per oscillation.

As this oscillatory motion scales down by $1/3$ each time and the initial
assumption was that $p_{0}$ corresponded to a “soft” leg, the scale of these
oscillations remains below that of hadronic sizes. Therefore the details of
the zig-zag string shape from this oscillatory behaviour are presumably not
significant and should not require detailed tracing for an accurate
hadronization picture. Thus, in this work we consider it sufficient to find an
average junction rest frame. Details on how the averaging treatment is carried
out are given in sec. 4.3. One should also note that though these diagrams
show a spatial representation of the strings for schematic simplicity and
visualisation purposes, in practice we handle everything from JRF calculations
to fragmentation in momentum space.

To construct the average, we will need the junction velocity during each
oscillation. This can be determined without an explicit boost to the Ariadne
frame and also generalised for massive endpoints. Consider a soft parton with
momentum $p^{\prime}_{0}$ in the Mercedes frame of a junction system (when
such a frame exists). Given the equation of motion
$\mathrm{d}p/\mathrm{d}t=-\kappa$, at time $p^{\prime}_{0}/\kappa$ the
endpoint will lose all its momentum and will change its direction of motion
back towards the junction, however in this frame the junction itself remains
at rest for now. This Mercedes frame will remain the JRF up until time
$2p^{\prime}_{0}/\kappa$, when the soft endpoint hits the at-rest junction,
with 3-momentum $-\vec{p\,}_{0}^{\prime}$ (irrespective of mass). This
momentum, $-\vec{p\,}^{\prime}_{0}$ will then dictate the direction of the
pull on the junction contributed by the string piece connected to the soft
parton. To calculate the pull from the other two junction legs we consider the
momenta of these legs at time $2p^{\prime}_{0}/\kappa$, by which time the
partons will have lost $2p^{\prime}_{0}$ momentum in their respective
directions of motion. Given these updated momenta for each junction leg, one
can construct a new Mercedes frame (assuming it exists), which then defines
the JRF for times after $2p^{\prime}_{0}/\kappa$ up till the next oscillation
begins. For massive oscillations, this will rapidly approach a pearl-on-a-
string scenario, which we discuss below.

####  3.2.2 Pearl-on-a-String

The next case to consider is a massive endpoint that is sufficiently soft that
a boost to the Mercedes frame does not exist. In such scenarios, the solution
to $f_{23}(|\vec{p\,}^{\prime}_{1}|)=0$ would return a negative-valued
$|\vec{p\,}^{\prime}_{1}|$ which of course is unphysical. Without resorting to
testing for negative solutions of $f_{23}(|\vec{p\,}^{\prime}_{1}|)=0$, such
configurations can be identified by looking at the rest frame of each massive
parton. Whilst in the rest frame of a given massive endpoint, if the opening-
angle between the other two legs is greater than 120°, a boost to a Mercedes
frame does not exist. One can easily convince oneself of this by considering
the boost required to reduce the angle of the other two legs to 120°, i.e.
$f_{23}(|\vec{p\,}^{\prime}_{1}|)=0$, and the resulting momenta of the massive
endpoint. Though not immediately obvious, this property is unique to at most
one parton in any given 3-parton configuration. In principle, as we assume
massless gluons these cases should only occur with endpoint partons. In such
cases, we expect the junction to get “stuck” to the massive endpoint, and
consequently one can map the junction motion by considering the massive parton
motion. When the junction becomes bound to the massive quark in this way, we
call that quark a “pearl-on-a-string”.

For the purpose of CR, simply using the initial rest frame of the pearl quark
is sufficient as at early times this will be the JRF, which has now been
explicitly implemented in Pythia 8.311 onwards. For fragmentation purposes we
wish to map out the junction motion. To do so, let us consider the simplest
case; a junction with a single massive endpoint and two hard massless
endpoints such that we have a pearl-on-a-string configuration. As with mapping
the oscillatory motion of a massless endpoint, here it is again instructive to
work with the equations of motion for this configuration in the Ariadne frame
with respect to the soft massive parton. Given the string tension $\kappa$,
the initial magnitude of the massive parton 3-momentum $p_{0}$, and the
position $x(t)$ and velocity $v(t)$ of the massive parton, we can write the
massive-parton momentum as a function of time,

$p(t)=p_{0}-2\kappa x(t)=\frac{mv(t)}{\sqrt{1-v(t)^{2}}}~{}.$ (14)

Note the factor of 2 here as the junction/pearl are connected to two other
string segments. Using the equations of motion of relativistic strings, eq.
(14) is derived by considering the energy gained by the string
$\mathrm{d}E_{s}$, over time $\mathrm{d}t$, given the massive parton moves
some distance $\mathrm{d}x$. Detailed derivations can be found in Appendix B.
This gives a maximal distance of $x_{max}=p_{0}/2\kappa$, at which we expect
the velocity to approach zero. Rearranging eq. (14) for the massive parton
velocity, we get the following:

$\frac{\mathrm{d}x}{\mathrm{d}t}=\frac{1}{\sqrt{1+\frac{m^{2}}{(p_{0}-2\kappa
x(t))^{2}}}}~{}.$ (15)

This differential equation is non-trivial to solve as the time dependence is
not straightforward, thus it is not very useful in this form for a practical
implementation. Instead of an exact solution, an approximation of the time-
dependent velocity can be used in its place. By numerically solving eq. (15),
the general shape of the solution can be approximated by a simple exponential,

$\frac{\mathrm{d}x}{\mathrm{d}t}=v_{0}\exp\left(\frac{-1.75t}{m}\right)~{},$
(16)

with $m$ being the mass of the pearl, $v_{0}$ the initial velocity in the
Ariadne frame, and using the assumption that $\kappa\sim 1$ GeV/fm. A
comparison of the solution $x(t)$ between the numerical solution and
approximate functional form is given in fig. 7, showing the motion of a pearl
quark over a time-scale of hadronization ($\sim$2 GeV) for several different
$p_{0}$ and $m$ values. Given this non-uniform motion of the junction bound to
a pearl with a decelerating velocity, the resulting string shape is curved and
appears as shown in the right image of fig. 8.

Figure 7: The pearl motion for an up/down, charm and bottom quark with
constituent masses, considering initial pearl velocities of $\frac{1}{2}$,
$\frac{1}{3}$ and $\frac{1}{10}$. Plots show the numerical solution to eq.
(15) (solid line), the simple approximation of the pearl motion given eq. (16)
(dashed line), and the motion of a gluon-kink on a $q\bar{q}$ string (dotted
line).

While this may seem like an oversimplification, we do believe it would be
overkill at this point to develop a more complicated functional form, also
just given the uncertainties, e.g., on the precise value of the string tension
to use, and the fact that the real physics is that of a curved string with
non-trivial time dependence which is not well understood. Rather than
attempting to fragment the actual curved string, we look to reuse techniques
for fragmenting dipole strings again in this context. The key behaviours we
wish to model in this pearl motion is the position (in momentum space) of the
junction baryon and some sense of the pearl $p_{\perp}$ that is propagated
along the strings to create the curve. The simple form of eq. (16) should be
sufficient for our purposes.

For the case of light quarks with sufficiently large $p_{0}$, one can see in
fig. 7 that the motion somewhat mimics that of a gluon kink on a $q\bar{q}$
string. In the Lund string model, gluons form form transverse excitations (or
kinks) on strings. These initially form a sharp corner on the string. Once
they have lost their momentum to the string, a flat so-called “region” forms
on the string, spanned between two kinks travelling in opposite directions.

A massless gluon with initial momentum $|\vec{p}_{g}|$ and initial velocity of
1 will lose its momentum after a time $|\vec{p}_{g}|/2\kappa$. Once it loses
all its momentum to the string, the flat region forms and the effective gluon
velocity becomes zero (there is no longer a force acting on the string in the
direction of the original gluon). The general shape of such a gluon-kink can
be seen in the left image of fig. 8. For a detailed description of string
regions and gluon kinks see Sjostrand:1984ic ; Sjostrand:1984iu . Given this
flat region formed by soft gluons, one naturally can make a comparison to the
curve in the pearl-on-a-string case; the pearl-on-a-string essentially acts as
a massive gluon.

Hence for light-quark cases with sufficiently large $p_{0}$, we may
approximate the massive pearl as a massless gluon, in order to mimic the
curved string and the $p_{\perp}$ propagated to the string from the pearl
quark. As such we can recycle the standard $q-g-\bar{q}$ fragmentation
procedure. In the following we will explain the procedure for the case of a
junction with quark endpoints; the same obviously holds for an antijunction
system with ($q\leftrightarrow\bar{q}$).

Figure 8: Left: a gluon kink on a $q\bar{q}$ string, forming a flat region
after the gluon momentum has been used up. Right: schematic of the junction
moving whilst “stuck” to a massive pearl, producing the shape of the pearl-on-
a-string scenario. Note this is by no means a to-scale diagram but merely to
show the general concept of the junction motion and the curvature of the
string. Each sketch shows the string configuration at two different times,
with black arrows indicating direction of motion. Both sketches are shown in
the Ariadne frame, the left-hand one with respect to the gluon, the right-hand
one with respect to the pearl.

In the gluon-pearl approximation, the massive quark and junction are replaced
by a massless gluon, $g_{q}$, with energy and momentum of $p_{0}$ in the
Ariadne frame. The excess energy from the mass of the pearl is “stored” at the
junction itself. From here the $q-g_{q}-q$ string can be mapped onto a
$q-g-\bar{q}$ string by taking $q\rightarrow\bar{q}$ for one of the junction
legs, and then standard dipole fragmentation procedures can be used. To
simplify the treatment, we initially force fragmentation from only the
$\bar{q}$ end of the $q-g_{q}-\bar{q}$ string. For each hadron produced from
this $\bar{q}$ string end, the conjugate is taken ($q\leftrightarrow\bar{q}$)
which maps it back to the initial $q-g_{q}-q$ system. This procedure is
continued until a string break steps over the junction, at which point the
pearl quark and the excess energy from the pearl quark mass is gained, forming
the junction baryon. As with standard junction fragmentation, the energy of
the produced hadrons is used as a measure of when the junction is reached.
Once the junction baryon has been made, all that remains is a typical
$q-\bar{q}$ string for which standard fragmentation is carried out including
the left-right randomization of string ends.

Notably this gluon-pearl approximation only holds for light-quark cases with
sufficiently large initial velocity $v_{0}$, such that the change in velocity
$\Delta v$, is approximately 1/2 over some reasonable hadronization time $\sim
2$ GeV. It is evident in fig. 7 that the approximation no longer holds for
heavier flavours or quarks with sufficiently small initial velocity. Let us
consider the limit where $p_{0}$ approaches zero and the limit of an
infinitely heavy quark. In such cases the pearl does not move and no
$p_{\perp}$ is propagated along the string. As such, we would expect the
standard junction hadronization method of fragmenting each leg towards the at
rest junction to work well. For large mass quarks or light quarks with low
non-zero $v_{0}$, although the junction motion is non-uniform, the change in
velocity is small. Thus in these cases one can use a perturbed JRF, which
simply takes the average velocity according to the approximate functional form
in eq. (16) over some hadronization time-scale. From here one can use the
standard junction fragmentation framework given this perturbed JRF.

There is of course ambiguity in which treatment to use in cases with
intermediate $p_{0}$ for light quarks and for charm quarks. Hence to create a
smooth transition between the two treatments, a probabilistic choice can be
made based on $\Delta v$ over hadronization times, with velocity changes of
zero and $1/2$ corresponding to the standard junction fragmentation and gluon-
pearl approximations respectively.

###  3.3 Gluon Kinks on Junction Motion

As with a simple $q\bar{q}$ string, each junction leg can also have any number
of intermediate gluons (kinks) between the junction and the endpoint. In this
section we will consider the effects of intermediate gluons on junction
motion. First we will look at a simplified model of only considering the pull
from the gluon extending out from the junction. This simplification is what is
implemented in Pythia 8.311 as outlined in sec. 4. Then we will allow the
gluon kinks to propagate back towards the junction, and the effect on the
junction motion as the kink passes through the junction. This is proposed
simply as a theoretical model and not yet implemented into Pythia as a general
description of multiple gluon kinks has not yet been solved. However we do not
expect the corrections due to gluon kinks propagating through the junction to
be large.

Figure 9: The example setup considered throughout sec. 3.3, with a single
gluon kink junction leg 3 and all massless parton. This shows the string
configuration at early times, depicted in the JRF of the first parton on each
leg, with the coloured dashed lines representing string segments.

We first consider the simplest example; a junction system with a single gluon
kink on one leg, and all partons massless. To examine the junction motion of
such a configuration, we consider the pull on the junction at different times.
As information is propagated along the string at the speed of light, the
initial junction velocity will be dictated by the first parton on each
junction leg, and the pull from partons further out along the leg would only
contribute to the junction motion at later times.

The string configuration of such a junction system with a single gluon kink in
its initial JRF is illustrated in fig. 9, where this early-time junction
motion is dictated by the partons labelled $q_{10}$, $q_{20}$, and $g_{30}$.
This frame is expected to remain the JRF until the first of these partons
loses its momentum, which is determined by the string equation
$dp/dt=-\kappa$. This equation of motion assumes a connection to a single
string with tension $\kappa$, however gluons are colour octets and therefore
connected to two other colour charges via string pieces, thus they lose
momentum twice as quickly compared to endpoints. The case of an endpoint being
the softest parton results in endpoint oscillations as described for massless
partons in sec. 3.2.1. Here let us instead consider the alternative scenario
where the gluon is the first to lose its momentum, i.e.,
$|\vec{p}_{g_{30}}|/2<|\vec{p}_{q_{10}}|$ and
$|\vec{p}_{g_{30}}|/2<|\vec{p}_{q_{20}}|$. The gluon will lose its momentum at
time $|\vec{p}_{g_{30}}|/2\kappa$, at which time $q_{10}$ and $q_{20}$ will
have lost momentum $|\vec{p}_{g_{30}}|/2$ in their respective directions of
motion.

Once the gluon has lost its momentum, the next parton on that leg, $q_{31}$,
will determine the pull on the junction, along with the reduced momenta of the
other two legs. Given the momentum of $q_{31}$ and the reduced momenta of
$q_{10}$ and $q_{20}$, one can construct a new JRF valid for times after
$|\vec{p}_{g_{30}}|/2\kappa$.

This procedure can be easily generalised to consider multiple gluons on each
junction leg by iteratively stepping outwards on each junction leg, and this
is the model we use for the implementation in Pythia which we elaborate on in
sec. 4.

Though this modelling of gluon kinks is adequate, similarly to expecting
endpoint oscillations one would expect the gluon kink to propagate back
towards the junction once it has lost its momentum, the behaviour of which is
ignored in the above description. Below we consider, for theoretical reference
without a full-fledged implementation, the effect on junction motion when
mapping out these gluon kinks in full detail.

Figure 10: For the junction configurations portrayed in fig. 9, each image
here shows which partons dictate the junction motion and the direction of the
pull from each such parton over a given time interval. From left to right, the
time intervals depicted are 0 to $|p_{g}|/2\kappa$, $|p_{g}|/2\kappa$ to
$|p_{g}|/\kappa$, and after $|p_{g}|/\kappa$. Note these momenta directions
are all shown relative to the initial JRF, and such arrows simply represent
the direction of the pull on the junction and not the magnitude.

In the same scenario depicted in fig. 9, at time $|\vec{p}_{g_{30}}|/2$ the
gluon $g_{30}$ will have lost its momentum, after which the gluon kink begins
to propagate back towards the junction. From here it will take another
$|\vec{p}_{g_{30}}|/2$ for the gluon kink to return and “hit” the junction.
This means that only after time $|\vec{p}_{g_{30}}|/\kappa$ will the junction
feel the effects of the kink propagating through the junction. Thus we can
split up the junction motion into three relevant time intervals relative to
the initial JRF; between 0 to $|\vec{p}_{g_{30}}|/2\kappa$,
$|\vec{p}_{g_{30}}|/2\kappa$ to $|\vec{p}_{g_{30}}|/\kappa$, and after
$|\vec{p}_{g_{30}}|/\kappa$.

Fig. 10 depicts each of these time intervals and the respective partons which
dictate the pull on the junction. As described prior, between times 0 to
$|\vec{p}_{g_{30}}|/2\kappa$ the JRF is defined by $q_{10}$, $q_{20}$, and
$g_{30}$ (left image of fig. 10), and at times $|\vec{p}_{g_{30}}|/2\kappa$ to
$|\vec{p}_{g_{30}}|/\kappa$ the JRF is defined by $q_{10}$, $q_{20}$, and
$q_{31}$ (middle image of fig. 10).

After time $|\vec{p}_{g_{30}}|/\kappa$, the gluon kink “hits” the junction and
propagates through to the other side with momentum $-\vec{p}_{g_{30}}/2$,
whilst the pull from parton $q_{31}$ remains the same. Thus in order to
determine the combined effect on the junction due to pulls from $q_{31}$ pull
and the gluon kink, one simply translates the momentum vector of $q_{31}$
along the direction of the propagating kink. This translation of the $q_{31}$
momentum is shown by the grey arrow in fig. 11.

Figure 11: Given the configuration set up in fig. 9, here shows a sketch
demonstrating the calculation of the pull leg 3 exhibits on the junction at
times after $|\vec{p}_{g}|/\kappa$.

As all partons in this example are assumed to be massless, the endpoints and
the kink move at the speed of light, and thus will expand out along the unit
circle. This means that the intersection of this kink-transported $q_{31}$
vector with the unit circle now can be used to define the overall direction of
pull on the junction from leg 3. This effective pull is labelled by
$q^{\prime}_{31}$ in fig. 11, which extends outwards from the centre of the
unit circle to the circle intersection with the kink-transported $q_{31}$.

Given an angle of $\theta$ between momenta of $g_{30}$ and $q_{31}$ in the
initial JRF, the angle between the initial $g_{30}$ direction and
$q^{\prime}_{31}$ will be $2\theta$. Hence at times after
$|\vec{p}_{g_{30}}|/\kappa$, the JRF is determined by momenta $q_{10}$,
$q_{20}$, and $q^{\prime}_{31}$ as shown in the right image of fig. 10. Note
here we have also assumed the endpoints are sufficiently hard such that there
are no endpoint oscillation effects.

Though the above provides a prescription for the treatment of a single soft
gluon on one junction leg, the calculation of these new pull vectors and their
associated times becomes increasingly complex when including variables such as
multiple gluons, massive endpoints, and/or soft endpoints. Hence we have not
developed this model to the point of being able to implement it into the JRF-
finding system as of Pythia 8.311. One would also expect that the other half
of the gluon momenta that propagates out towards the endpoint of the junction
leg would eventually return towards the junction. However due to hadronization
timescales, the junction motion at late times are not expected to contribute
significantly if at all, and such detailed modelling may well be overkill.

##  4 Implementation

In the preceding section we explored how to map the junction motion out over
time, including oscillation effects, pearl-on-a-string scenarios, and
junctions with gluon kinks. Here we explore how to apply these ideas to a
practical hadronization model.

For pearl-on-a-string cases (when a junction is effectively bound to a
soft/slow endpoint), the implementation follows the procedure described in
sec. 3.2.2. The approximation of using a gluon to represent the soft/slow
endpoint turns out to only be practically useful in a somewhat limited set of
circumstances; we only use it in cases when gluon-kinks do not cause changes
in the junction motion over a time scale characteristic of hadronization (see
below). This means that only the first parton on each junction leg dictates
the junction motion, whether it be because the non-pearl junction legs are
sufficiently hard or are endpoints. We have not generalised the gluon
approximation to junction topologies with soft gluon kinks that cause changes
in the junction motion as in such cases there is no longer a straightforward
method to approximate the massive pearl as a massless gluon. The gluon
approximation is also only considered if the pearl has at least its
constituent mass in order to protect against incorrectly assigned quark
masses.

The scale that determines the softest kink for which the gluon-approximation
may be used is set by the parameter StringFragmentation:pNormJunction, which
we label $p_{norm}$ below. This should reflect an average (inverse) time for
string breaks to occur (i.e.
$\langle\Gamma\rangle=\langle\kappa\tau\rangle\sim 2$ GeV).

When standard junction fragmentation is used instead (the majority of cases),
the JRF is needed to create the fictitious legs; this is especially important
to predict the junction baryon kinematics. However, to deal with junction
systems with gluon-kinks where the junction motion changes over time, instead
of mapping out the changing junction motion in detail (and e.g. dividing it
into different string worldsheet regions), we here simply use a single average
JRF. This notion of an average JRF was also used in the previous
implementation Sjostrand:2002ip , however here we provide a new method to
determine the average junction motion in way that is more stable and reliable
when handling gluon kinks and soft/slow endpoints. In the following
subsections we outline the new iterative procedure to calculate the average
JRF, and explain the prescription used for fragmenting such systems.

###  4.1 Time dependence of Junction Motion

We first briefly outline the general idea of the procedure implemented in
Pythia 8.310 and prior (with full details in Ref. Sjostrand:2002ip ) — the
“old” procedure — before turning to the new one.

In the old JRF finding procedure, weighted sums of the four-momenta of the
partons on each junction leg are made, which are called “pull vectors”. The
JRF is defined as the frame in which these pull vectors form a Mercedes
configuration. The weightings used in the construction of these pull vectors
depend on parton energies and are thus frame-dependent. The determination of
the JRF is therefore an iterative procedure. For a given set of pull vectors,
one boosts to the Mercedes frame defined by those pull vectors; one then
updates the pull vectors. These updated vectors are not necessarily in a
Mercedes configuration; one then moves to the Mercedes frame defined by the
updated vectors, and iterate until, ideally, the procedure converges.
Convergence is not guaranteed however, especially if/when large-mass pull
vectors result due to the summation of four-momenta. Indeed, the iterative
procedure fails in around 10% of minimum-bias events at LHC energies, in which
case the procedure reverts to the centre-of-mass frame as a fall-back frame in
which to the fragmentation instead. Moreover, the procedure assumes the JRF
must be of Mercedes type, and allows root-finding for the Mercedes frame to
return an unphysical answer in the would-be pearl-on-a-string cases (when the
Mercedes frame does not exist). This prescriptions was largely formulated with
high-energy legs in mind, and does work well in that context however falls
short with many soft gluon kinks and soft endpoints.

The procedure we propose here has been formulated to be valid for arbitrary
combinations of gluon kinks and endpoint masses. It does not presuppose the
existence of a Mercedes frame nor does it rely on convergence of an iteration.
Instead, it finds a time-ordered sequence of well-defined JRFs each of which
is valid in a given time window, and then makes a time-weighted average over
these successive JRFs to find an overall average JRF. To determine the motion
of the junction, the procedure steps sequentially through the partons on each
junction leg. As per the simple description of gluon kinks in sec. 3.3, at
early times the partons immediately nearest to the junction will dominate the
pull on the junction. After the momentum of the nearest parton on a leg has
been depleted, the next parton on that junction leg will take over in
dictating the junction motion, and so on. To carry out this procedure, we need
to keep track of several pieces of information per iteration; the junction
velocity, the time interval the velocity is relevant for, and the four-vectors
that will dictate the next iterative JRF.

To simplify the language used below, we freely use momentum magnitude as a
measure of time. This is justified as the time it takes each parton to lose
its momentum will be proportional to the 3-momentum magnitude, $|\vec{p\,}|$,
of the parton according to string equation of motion $dp/dt=-\kappa$, with an
additional factor of 2 for gluons. We also simplify notation by distinguishing
between endpoint parton momenta (inclusive of (anti)quarks and (anti)diquarks)
and gluon momenta with subscripts $q$ and $g$ respectively.

Junction equations of motion and the calculation of junction velocities were
discussed in sec. 3. For a given set of four-vectors, one finds either the
Mercedes frame should it exist, or one uses the pearl-on-a-string notion and
find the average velocity using the approximation in eq. (16). We expect each
JRF to remain the rest frame until either a gluon depletes its momentum, or an
endpoint parton oscillates and returns to hit the junction. For a given
configuration, the time until the next change in junction velocity is given by
whichever has the smallest momentum among the three partons that are currently
adjacent to the junction. This “currently smallest momentum” is labeled
$p_{\rm small}$.

For a Mercedes-frame topology, $p_{\rm small}$ is defined as the smallest of
any $|\vec{p}_{g}|/2$ or $2|\vec{p}_{q}|$ in the given three-parton
configuration. The factor of $1/2$ on $|\vec{p}_{g}|$ is to account for a
gluon being connected to two string pieces, hence it loses energy twice as
fast and half of its momentum propagates outward, away from the junction,
while the factor of 2 on $|\vec{p}_{q}|$ helps to account for endpoint
oscillations as will be described below. It also accounts for the fact that
after an endpoint has lost all of its energy, it takes the same amount of time
again for that information to propagate back to hit the junction. For pearl-
on-a-string cases, we define $p_{\rm small}$ as the smallest momentum of the
two non-pearl junction legs in the Ariadne frame. As this is with respect to
the Ariadne frame and not the JRF, this time is multiplied by a
$\gamma-$factor to transform from the Ariadne frame to the perturbed JRF. As
these $p_{\rm small}$ times (in both the Mercedes and pearl-on-a-string cases)
are measured in each iterative JRF, an additional $\gamma$-factor is used to
translate these associated times back to the laboratory frame.

Once the junction velocity (in the lab frame) and the associated time is
calculated for a given set of four-vectors, the momenta are then updated for
the next time interval. For pearls, the momentum at time $p_{\rm small}$ is
simply determined by the velocity at time $p_{\rm small}$ according to eq.
(16). Otherwise using the momentum-loss relation
$\mathrm{d}p/\mathrm{d}t=-\kappa$ again, each parton will lose momentum at a
constant rate irrespective of the parton mass. This means after a time $p_{\rm
small}$, the partons will now have updated 3-momenta of
$\vec{p}\,^{\prime}_{i}=\hat{\vec{p}}_{i}(|\vec{p}_{i}|-p_{\rm small})$, with
their energies scaled accordingly to preserve their mass. If the parton that
defines $p_{\rm small}$ is an endpoint, this momentum-updating mechanism
naturally incorporates oscillations about the junction as it will update the
3-momentum to be $-\vec{p}_{q}$. If the parton that defines $p_{\rm small}$ is
a gluon, one simply updates the pull from that junction leg by stepping to the
next parton on said leg.

Putting the above together, we can now construct a full sequential procedure
for finding the junction motion over time. The steps of this process are as
follows, with the initial 3-parton configuration being the first parton on
each junction leg. A schematic example of a corresponding sequence is given in
fig. 12.

  1. 1.

Check whether we have a pearl-on-a-string configuration or a Mercedes type
JRF. Compute and store the velocity $v_{i}$ with respect to the initial frame
of reference, and boost to this frame.

  2. 2.

Find $p_{{\rm small}_{i}}$ for the given 3-parton configuration.

  3. 3.

Update the four-vectors that dictate the junction motion, by stepping to the
next parton on the leg associated with $p_{{\rm small}_{i}}$ if possible, else
update momenta according to the time $p_{{\rm small}_{i}}$.

  4. 4.

Boost the system back to the initial frame and repeat the process till either
the sum of all $p_{\rm small}$ exceeds a threshold value of $p_{\rm max}$
(definition proceeds in the next section), or till two endpoints have been
reached.

Figure 12: Illustration of the JRF averaging procedure, by stepping through
individual JRF frames sequentially. The three-momenta of partons along each
leg are lined up, stepping outwards along the junction legs. The black arrows
highlight the momenta that dictate junction motion. The left image shows the
first step in the sequential procedure: finding the initial JRF, dictated by
the first parton on each leg. The middle image illustrates the combination of
steps 2 and 3 of the sequential procedure procedure. In this example $g_{30}$
defines $p_{\rm small}$. After that amount has been lost from each leg,
highlighted in cyan, the pull for leg 3 is now determined by the next parton
on that leg, $g_{31}$, and the momenta remaining on the other legs are
updated. The right image depicts the start of the next iteration, showing the
next JRF given the updated momenta for each junction leg. Note that angles and
vector lengths here are not to scale and are only intended to give a visual
representation of the sequential procedure.

###  4.2 Handling of Exceptional Topologies

Though the above procedure is generally stable, there are still a few
topologies we have made exceptions for: massive gluons and collinear massless
partons.

Though the default gluon mass in Pythia is zero, non-zero mass gluons can be
encountered at times either from user input or the clustering of two or more
nearby gluons into a single small-mass gluon. Should the Mercedes frame not
exist for a configuration due to a massive gluon being soft, these are not
treated as pearl-on-a-string cases. Instead the contribution from the soft
gluon to the JRF is considered negligible and we step to the next parton on
that leg.

The other special case considered is collinear massless partons, which is
mostly expected to occur with massless gluons. In such cases a boost to a
Mercedes frame will never be possible, nor does it make sense to think about a
pearl-on-a-string. Perfectly collinear partons are not expected to be
encountered often in Pythia event generation, however numerical precision of
nearly collinear partons can lead to issues in finding a Mercedes frame given
the root-finding procedure used. Additionally we wish to ensure stability if
given unphysical user-inputted parton configurations. For a pair of exactly
collinear partons in a three-parton configuration, the junction motion becomes
ill defined however one would expect the junction to be highly boosted in the
direction of the collinear pair. Hence the following fail-safe is implemented
not so much to describe the exact physics of the scenario, but to provide some
good approximation to reflect the boost in the collinear pair direction and to
give stability to the procedure.

To handle these configurations, we form a four-vector in the direction of
motion of the collinear pair and assign it a diquark mass, and then use this
diquark-type velocity as the junction velocity. To construct this, in the
centre-of-mass frame of the 3-parton configuration, the summed momenta of the
collinear partons defines the direction and energy of this diquark-mass four-
vector. The three-momentum magnitude is then fixed by the constituent mass of
diquark $ud_{0}$ which is the lightest diquark according to the particle data
table in Pythia. From here the velocity of this four-vector defines the JRF,
and the three-momentum magnitude defines $p_{\rm small}$. The momenta used in
the iterative JRF are then updated to the next parton on all junction legs.

If any of the partons were already endpoints and hence one cannot step out
along the leg further, we simply stop the iterative procedure here. Should the
collinear pair have insufficient energy to form a diquark mass and both
partons are gluons, we simply step to the next parton on each of these legs
and continue with the next iteration. If there is insufficient energy and one
of the collinear partons is an endpoint, no junction velocity is recorded for
this iteration and the iterative procedure is stopped here. In the rare case
we encounter this scenario and there were no previous iterative JRFs found,
i.e. the soft collinear partons are the first partons on their respective
junction legs and at least one is an endpoint, then we resort to a fail-safe
of defining the JRF as the centre-of-mass frame. Indeed junctions should not
be forming directly between collinear partons, however nonetheless we have
ensured to protect against such occurrences and to ensure procedural
stability.

###  4.3 Average Junction Rest Frame

Using the set of sequentially calculated velocities $\vec{v}_{i}$, according
to the above procedure, an average JRF, and hence the associated average
junction velocity, $\vec{v}_{\rm jun}$, can be determined. In the following,
the first iteration of the sequential procedure is marked by $i=1$. Since
fragmentation will gradually happen, we expect that the pull on the junction
at early times will be more important than those at late times to determine
its motion over the time scales relevant for hadronization. We introduce a
time-dependent exponential decay to weight each $\vec{v}_{i}$. As explained
above, the length of time each JRF is relevant for is given by $p_{\rm
small}$. Importantly this time measure is defined in the JRF, and thus must be
translated to the initial frame of reference by a Lorentz factor,
$\gamma_{i}=1/\sqrt{1-|\vec{v}_{i}|^{2}}$. This allows the time interval for
junction motion $\vec{v}_{i}$ to be defined from times $p_{i-1}$ to $p_{i}$,
where $p_{i}$ is a sum of times defined by
$p_{i}=\sum_{j=0}^{i}\gamma_{j}p_{{\rm small}_{j}}$ with $p_{0}=0$. Using
this, the calculation of the average junction velocity is given by,

$\vec{v}_{\rm jun}=\frac{\sum_{i=1}^{i_{\rm
max}}\vec{v}_{i}(e^{-p_{i-1}/p^{\prime}_{\rm
norm}}-e^{-p_{i}/p^{\prime}_{norm}})}{1-e^{-p_{i_{\rm max}}/p^{\prime}_{\rm
norm}}}.~{},$ (17)

where the parameter $p\textquoteright_{\rm norm}$ acts as a reference scale
characteristic of the time scale of the hadronization process. It is defined
via the parameter $p_{\rm norm}$ (set by StringFragmentation:
pNormJunction), which is assigned the default value of 2 GeV. If we want
$p_{\rm norm}$ to define a (proper) time in the JRF, we need to redefine this
normalisation parameter for the initial frame of reference in which we are
calculating the average junction velocity, which we here we call $p_{\rm
norm}^{\prime}$. To do so, we consider the sum of times $p_{{\rm small}_{i}}$
up till they add to $p_{\rm norm}$, then incorporate the changing junction
motion by a Lorentz factor. The summation of these $\gamma-$scaled $p_{\rm
small}$ values then defines $p^{\prime}_{\rm norm}$,

$p^{\prime}_{\rm norm}=\sum_{i=1}^{N}\gamma_{i}p_{{\rm
small}_{i}}+\gamma_{N+1}(p_{\rm norm}-\sum_{i=1}^{N}p_{{\rm small}_{i}}.$ (18)

The second term in eq. (18) ensures the sum of $p_{\rm small}$ add to $p_{\rm
norm}$ exactly, where $N$ is defined such that $\sum_{i=1}^{N+1}p_{{\rm
small}_{i}}>p_{\rm norm}$ and $\sum_{i=1}^{N}p_{{\rm small}_{i}}\leq p_{\rm
norm}$. In the case $p_{{\rm small}_{1}}>p_{\rm norm}$, we simply have
$p^{\prime}_{\rm norm}=\gamma_{1}p_{\rm norm}$. We also use the parameter
$p_{\rm norm}$ to control the value of $p_{\rm max}$ which dictates the
stopping point in the sequential JRF averaging procedure. We choose $p_{\rm
max}=5p_{\rm norm}$ so that we keep calculating junction velocities beyond the
average hadronization time whilst having these JRFs heavily suppressed by the
exponential decay in eq. (17).

Once the average JRF has been constructed, the fictitious endpoints are formed
by summing the momenta on each junction leg and reflecting this summed momenta
on the other side of the junction, after which standard fragmentation can then
be carried out, as described in Sjostrand:2002ip ; Bierlich:2022pfr .

##  5 Results

We begin this section by examining the theoretical expectations from the
updated modelling of junctions. We then inspect the effects in hadron-
collision events by comparing to both experimental data and to the previous
junction modelling. In the following, the treatment described in this paper
(and implemented in Pythia 8.311) is labelled “new” while the treatment in
Pythia 8.310 is labelled “old” Sjostrand:2002ip ; Christiansen:2015yqa .

###  5.1 Theoretical Implications

We have made three key alterations to the junction modelling implementation in
Pythia 8.311; the $\lambda$-measure used in the QCD CR algorithm, the junction
rest-frame finding procedure, and the use of a gluon approximation for
fragmentation of (a subset of) junction systems with soft light-flavour legs.

String-length measure changes: We first verify that the change to the new
$\lambda$-measure, eqs. (4) and (7), does not cause unexpected large effects
in comparison to the previous form eq. (3). In fig. 13 we show the ratio of
all junction baryons to the total baryon yield (solid lines), as well as
specifically the heavy-junction baryon to total heavy-baryon yield (dashed
lines), as a function of baryon $p_{\perp}$ (left) and rapidity (right).

Figure 13: Junction-baryon to total baryon yield ratios, considering both
ratios of all flavours of baryons (solid lines) and heavy-flavour baryons only
(dashed lines). We compare both the old $\lambda$-measure for both the new and
old junction modelling, to the updated $\lambda$-measure in conjunction with
the new junction model. Simulations were for NSD events at $\sqrt{s}=7$ TeV,
looking at hadrons produced directly from hadronisation (i.e., with status
codes between 80 and 90).

This demonstrates that there are no major changes at the level of total
yields. It also highlights the importance of rigorous junction-motion
modelling as up to 50% of light-flavour baryons and up to 75% of heavy-flavour
ones are coming from junctions. We return to this below.

We note that we do expect small differences, e.g., due to the parameter
$j_{0}$ in eq. (7)222Set by ColourReconnection:junctionCorrection, with
default value 1.2., which modifies the probability of junction reconnections.
Although the same parameter is used in both the new and old
$\lambda$-measures, the impact of the parameter is not identical between them.

For completeness, although our study is focused on junction topologies we note
that changes to the string-length measure will also affect dipole connections,
which may in turn affect, e.g., heavy-flavour meson production and potentially
the frequency of usage of the ministring fragmentation Norrbin:2000zc
procedure. To facilitate investigations of differences, both forms of the
$\lambda$-measure are available in Pythia 8.311, controlled by the parameter
ColourReconnection:lambdaForm. For the remainder of the paper, when using the
new junction modelling in Pythia 8.311, we use the new string-length measure
(ColourReconnection:lambdaForm = 0 as of Pythia 8.311).

Junction Rest-Frame Determination: The changes to the JRF-finding procedure
comprise the new junction velocity-averaging method as well as taking into
account the pearl-on-a-string motion for soft legs. The former is relevant for
junction topologies with one or more gluons along the legs, and is difficult
to illustrate in a simple/clean way. We therefore defer discussion of these
cases to sec. 5.2 on experimental comparisons. Here, we focus on the extent to
which the JRF changes for pearl-on-a-string cases affects total hadron
multiplicities as well as the $p_{\perp}$ spectra of the junction baryon and
of the hadrons from the other junction legs.

Figure 14: A simple setup for a junction system with two massless junction
legs aligned along the z-axis and defined by $p_{2}$ and $p_{3}$, and a single
masive quark defined by $p_{d/b}$ with momentum along the positive x-axis. In
this setup, the massive quark is either a down or bottom quark with
constituent masses 0.33 GeV and 4.8 GeV respectively, and down quarks are used
for the massless quarks defined by $p_{2}$ and $p_{3}$. This configuration is
constructed in the Ariadne frame with respect to the massive quark, with
$p_{\perp}$ measured relative to the z-axis.

Consider the simple generic test case illustrated in fig. 14: a three-parton
junction configuration with two energetic massless quarks and a single massive
one, whose properties we can vary. The momenta of the partons are set up in
the Ariadne frame of the single massive quark. That is, the two massless
quarks are aligned to be back-to-back along the $z$-axis, and the momentum of
the massive quark is aligned with the positive $x$-axis. The energies of the
massless quarks are arbitrarily both set to 50 GeV, mainly just to ensure that
the massive quark is always the softest leg and that there are no soft effects
associated with the other two legs.

This particular configuration is chosen for convenience; given that the two
energetic legs are massless and back-to-back, a boost by $v=1/2$ will reduce
their opening angle to 120° , as required in the Mercedes frame. Thus, pearl-
on-a-string cases are easily identifiable and correspond to massive-quark
velocities of less than 1/2.

To further simplify these tests and isolate junction-baryon production more
cleanly, we also turn off diquark production from standard string breaks
(StringFlav:probQQtoQ = 0) and the popcorn mechanism which allows for meson
production from diquark endpoints (StringFlav:popcornRate = 0).

We consider two limiting cases for the massive quark; a light down quark with
a constituent mass of 0.33 GeV and a heavy bottom quark with a constituent
mass of 4.8 GeV.

First, let us discuss the expectations from the old JRF-finding procedure.
When the Mercedes frame does not exist (i.e., for massive-quark velocities
less than 1/2 in the Ariadne frame, a.k.a. pearl-on-a-string cases), the old
algorithm accidentally converges on the rest frame of the pearl quark. Hence,
for such cases the junction velocity (in the Ariadne frame) will simply be the
quark velocity. For Ariadne-frame quark velocities greater than 1/2, the
Mercedes frame exists and the old algorithm takes this as the JRF;
oscillations of soft legs around the junction are not mapped out and therefore
the junction velocity (in the Ariadne frame) will simply be 1/2.

The new modelling seeks to take into consideration the deceleration of the
pearl quark due to the connection to two other string pieces, determining an
average JRF which, for pearl-on-a-string scenarios, should have a lower
velocity (in the Ariadne frame) than that of the massive quark. For cases
where the Mercedes frame does exist, the junction velocity may also still be
lower in the new treatment, if the endpoint energy is low enough that
oscillations about the junction become relevant (taken into account in the new
treatment, ignored in the old one).

The old and new treatments will coincide when the massive-quark leg velocity
$v>1/2$ and it is hard enough that no oscillations need be considered. The
scale of “hardness” here when talking about oscillations is relative to
$p_{\rm norm}$ (StringFragmentation:pNormJunction), where the junction motion
is most heavily weighted for times within $p_{\rm norm}$. For the two limiting
cases we consider here, a $b$ quark with mass 4.8 GeV and a light quark with
mass 0.33 GeV, the no-oscillation limit (for the default value of $p_{\rm
norm}=2\,\mathrm{GeV}$) corresponds to quark velocities above 0.64 and 0.98,
respectively.

Hence, for the light-quark cases in figs. 15, 16 and 17, which are all
distributions as a function of the Ariadne-frame massive-quark velocity, the
convergence is not quite visible as it happens at the right edge of the plots.
We also note that all distributions in figs. 15, 16 and 17 are for the hadron
spectra given the massive quark leg does not fragment, or in other words the
massive-quark endpoint becomes contained within the junction baryon itself.

Figure 15: The average transverse momentum (top row) and multiplicity (bottom
row) of hadrons produced from the fragmentation of the two junction legs
defined by the massless quarks given the Ariadne frame setup described in fig.
14. These distributions are given as a function of the massive quark velocity
(in the Ariadne frame), where the massive quark is a down quark in the left
column of plots, and a bottom quark in the right column.

We will begin by making a comparison of the JRF finding procedures given the
standard junction fragmentation method (with the use of fictitious legs), then
proceed to compare to the gluon-approximation approach. For this test scenario
as per fig. 14, let us first let us examine the effect on the fragmentation of
the two more energetic junction legs defined by $p_{2}$ and $p_{3}$.

Since the boost between the Ariadne frame and the JRF is lower in the new
treatment than in the old (in which the massive-quark rest frame was used for
the JRF), the magnitudes of $p_{2}$ and $p_{3}$ will be less in the new
treatment. Thus when constructing the fictitious endpoints used to approximate
each leg as a dipole, the string lengths would be a bit smaller in the new
average JRF compared to if the rest frame of the massive quark were used.
Hence one would expect the new JRF-finding procedure to result in an overall
smaller number of hadrons produced from the two harder junction legs. This
result is evident in the bottom row of plots in fig. 15, which shows the
average hadron multiplicities from junction legs defined by $p_{2}$ and
$p_{3}$ as in fig. 14.

Not only would the lower junction velocity affect the hadron multiplicities,
it should also impact the $p_{\perp}$ distributions of the hadrons. The
construction of the dipole using the fictitious endpoint does not just
determine the string length but also dictates the string axis which defines
both the longitudinal direction of string breaks (according to the symmetric
fragmentation function in eq. (6)), and the axis relative to which $p_{\perp}$
kicks from string breaks are determined (according to the Schwinger mechanism
in (5)). The lower average junction velocity in the new treatment will result
in the fictitious dipoles of the massless quark legs being less boosted, hence
should result in a softer average $p_{\perp}$ distribution for the produced
hadrons. This effect can be seen in the top row of plots in fig. 15.

Figure 16: The average $p_{\perp}$ distributions for the junction baryon (top
row) and junction diquark (bottom row), given the setup in fig. 14, such that
the massive quark is contained within the junction baryon. The left and right
columns are for massive quark flavour of down and bottom respectively.

Figure 17: The distributions here are for the setup given in fig. 14 with a
down quark of mass 0.33 GeV for the massive quark. As a function of the
massive quark velocity the left panel shows the average $v_{x}$ of the
junction baryon produced, and the right panel shows the mass of the junction
diquark produced by the standard junction fragmentation procedure. The
junction baryons and diquarks considered here are only those which contain the
massive down quark in the initial setup.

Most interesting is of course the effect on the junction baryon itself.
Predicting this in the updated treatment is less straightforward however.
Instinctively one would expect a lower junction velocity to correspond to a
lower junction baryon velocity. However this effect is skewed by the old
treatment overpredicting the $p_{\perp}$ distribution of the more energetic
junction legs. When carrying out the standard junction fragmentation
procedure, after the two softest legs are fragmented, the momentum of the
junction diquark (see fig. 5 for reference) is determined using energy-
momentum conservation according to eq. (8). As the $p_{\perp}$ of the produced
hadrons is overpredicted when using the massive-quark rest frame, this is
compensated for by a softer junction diquark. Hence when looking at the
$p_{\perp}$ distributions of the junction baryon when the massive quark is a
bottom quark (top right panel of fig. 16), the new junction modelling actually
results in a slightly harder baryon despite having a lower junction velocity.

Interestingly, for the case of the massive quark being a light-flavoured quark
(with a constituent mass), we observe the opposite behaviour. However though
the $p_{\perp}$ of the junction baryon is greater as seen in fig. 16, the
direction of motion of the baryon is oppositely oriented to the junction
velocity itself. As evident in the left-hand panel of fig 17, the old
modelling results in the junction baryon traveling in the negative
$x$-direction, despite the JRF per the old modelling being the pearl-quark
rest frame, which in this example is moving in the positive $x$-direction.
Naturally, we would not expect the junction baryon to move in the opposite
direction to the junction itself; this unwanted behaviour in the modelling is
a consequence of the incorrect JRF being used. Using the light-quark rest
frame as the JRF overestimates the junction velocity, thus when fragmenting
the second-softest junction leg in this frame, the hadrons produced take too
much positive $p_{x}$, which then is compensated for when constructing the
junction diquark, and then the junction baryon.

Another contributor to this discrepancy is the negative-mass junction diquarks
in the old treatment, cf. the right-hand panel of fig. 17. This issue arises
from only using the energy as a measure of when to stop fragmenting the two
soft legs in towards the junction, and is avoided by the additional constraint
added in the new procedure that checks both the remaining energy of each
junction leg and the mass of the potential junction diquark. We note that the
positivity of the diquark masses in the new treatment also noticeably reduces
the number of errors encountered during the fragmentation of the last string
piece.

Overall we conclude that, for the standard junction fragmentation procedure,
using an average JRF instead of simply the early-time JRF (as per the old
modelling) results in lower hadron multiplicities and lower average
$p_{\perp}$ distributions (excluding the junction baryon). The effect on the
junction baryon itself is less straightforward but the treatment of the
junction diquark in the updated modelling is more stable and results in fewer
errors333Typically of the kind ‘‘StringFragmentation::fragment: stuck in
joining’’, which in Pythia 8.311 is now encountered approximately a factor 15
less (for dijet events at LHC energies)..

Figure 18: The $\Lambda/K^{0}_{S}$ ratio with respect to $p_{\perp}$ (left
panel) and $|y|$ (right panel), with data from CMS CMS:2011jlm . The models
shown here are using Pythia 8.311 with and without the gluon approximation
allowed for pearl-on-a-string cases. All events are 7 TeV NSD events, where
Pythia simulations have a lifetime cut of $\tau_{max}=10$ mm/c and no
$p_{\perp}$ cuts on final state particles.

Gluon Approximation for Light Pearls: Next we consider the effect of using the
gluon-approximation approach for fragmenting light-flavoured pearl-on-a-string
cases. Note here we neglect the bottom-quark case as the gluon approximation
only makes sense in light quark cases, cf. the arguments given in sec. 3.2.2.
The aim of the gluon approximation is to capture the effects of the motion of
the light pearl quark, which as seen in fig. 7 traverses a similar distance
over time pNormJunction as a massless gluon kink with the same 3-momentum.
This should ideally mimic the $p_{\perp}$ from the pearl being propagated out
along the two non-pearl junction legs, which in turn should mimic the curve in
the strings better than approximating each junction leg as half of a dipole
string.

As the $p_{\perp}$ in this case would be explicitly propagated along the
string by the gluon kink, one would expect a harder $p_{\perp}$ spectrum
overall in comparison to the standard fragmentation approach, which is evident
in the top left panel of fig. 15. The $p_{\perp}$ spectrum given the gluon
approximation still remains below that of the old modelling. Interestingly the
hadron multiplicity is higher than the standard fragmentation approach and
even becomes higher than the old modelling as seen in the bottom left panel of
fig. 15. Given a higher $p_{\perp}$ and $n_{\text{had}}$, this means reduced
momenta longitudinally. This is somewhat surprising and appears to be a
consequence of the modelling of gluon kinks. Whether this achieves the desired
result is somewhat difficult to determine as there is ambiguity in how well
one can approximate a massive pearl as a massless gluon. We retain this case
mostly out of theoretical interest and to allow to explore modelling
ambiguities.

As the gluon approximation aims to distribute the $p_{\perp}$ of the pearl
along the string, one may assume that the junction baryon itself would retain
less $p_{\perp}$. However interestingly the junction baryon $p_{\perp}$ is
predominantly harder in these cases compared to the standard fragmentation
procedure. From looking at the $p_{\perp}$ of the junction baryon alone, one
may initially assume the baryon becomes harder in the junction direction of
motion (the $x$-direction, in our example). However by looking at the velocity
in the $x$-direction of motion in the left panel of fig. 17, one can see that
it is actually reduced when using the gluon approximation. This means that the
$\langle p_{\perp}\rangle$ from string breaks is greater, however less biased
in the positive $x$-direction.

One would reasonably expect the standard procedure to bias the junction
direction of motion more in the positive $x$-direction, as the junction baryon
necessarily retains the momentum of the massive quark. In contrast, in the
gluon approximation, the $p_{x}$ of the pearl is explicitly propagated
outwards along the hard-leg string pieces, resulting in a junction baryon with
less positive $p_{x}$. It is therefore somewhat surprising that the
$p_{\perp}$ of the junction baryon in the gluon-approximation case is larger,
however this again is a result of the treatment of gluon-kinks and the
$p_{\perp}$ from string breaks. Comparatively, the old junction model produced
junction baryons predominantly moving in the negative $x$-direction, which is
opposite to the junction motion. This negative average $v_{x}$ is another
consequence of negative-mass junction diquarks in the old treatment.

Note the in the current implementation, the gluon approximation is only used
in cases where the pearl forms at initial time zero, meaning that only cases
with quark velocities of less that 1/2 in this Ariadne frame setup are
considered for the gluon-approximation approach. Though the treatment could in
principle be extended to pearls that form at later times which in turn would
smooth out these distributions, this is not considered in the current
implementation. This is partially as there is no clear approach to what would
be a reasonable gluon momentum to use in these cases, but also these gluon
approximation cases are not encountered often in real events which will be
further elaborated on below.

The strict cutoff for the gluon approximation at $v=1/2$ results in the curves
in fig. 15 and fig. 16 showing a sharp change at down-quark velocity of 1/2.
One should also keep in mind that the approximation of the massive pearl quark
as a massless gluon kink is not perfect. Though approximating the pearl as a
gluon results in an overall similar displacement of the pearl as shown in fig.
7, the velocities of a pearl quark and a gluon kink are considerably
different, with the pearl having a maximum initial velocity of 1/2 whereas a
massless gluon will have a velocity of 1. Given this discrepancy between a
realistic pearl quark and a gluon kink, it is somewhat difficult to tell
whether the gluon approximation or simply using an average JRF is a more
faithful approach in such cases.

Another factor to consider is that the gluon approximation is also only well
defined in a simple three-parton configuration as described above. In the case
of soft gluon kinks, there is no longer a distinct Ariadne frame to perform
the gluon construction in, nor would the pearl motion necessarily be well
mimicked by a gluon in such cases anyway. Given this constraint, when looking
at a hadron collision event such as a $\sqrt{s}=7$ TeV $pp$ collision, only
around 0.01% of junction baryons are formed using the gluon approximation. We
therefore expect this modelling to have minimal practical consequences.

Additionally, any observable effects of this modelling would be further
diluted by the portion of light-flavour baryons produced via standard string
breaks via diquark-antidiquark pair creation. Thus though this model provides
an in-principle elegant parallel between the pearl-on-a-string junction cases
and gluon kinks, in practice it is of quite limited use as these cases rarely
occur. This is evident in fig. 18 where there is no significant observable
impact on the $p_{\perp}$ distribution of either the $\Lambda/K$ ratio or the
$\Xi/\Lambda$ ratio with the inclusion of the gluon approximation. Hence, for
the comparisons in sec. 5.2, we do not include the gluon approximation.

###  5.2 Experimental Comparisons

Here we examine the results of the revised junction fragmentation modelling in
comparison to experimental data, with a particular focus on baryon-to-meson
ratios. Note that in the following we do not attempt a full retuning but leave
most parameters untouched relative to the “CR Mode 2” tune of the QCD CR model
Christiansen:2015yqa . This allows for a direct comparison between the old and
new junction modelling. We also will compare results to the Monash (2013) tune
Skands:2014pea , which provides a baseline model without junctions. A full
detailed tuning effort will be left for a future study.

Figure 19: The charged multiplicity average as a function of pseudorapidity on
the left CMS:2010tjh , and the right plot shows the average $p_{\perp}$ as a
function of charged multiplicity CMS:2010qvf . All events are 7 TeV NSD events
from CMS, where Pythia simulations have a lifetime cut of $\tau_{max}=10$ mm/c
and no $p_{\perp}$ cuts on final state particles.

Only two of the CR Mode 2 baseline parameters were altered:

  * •

MultiPartonInteractions:pT0Ref was increased slightly, from 2.15 GeV to 2.25
GeV, to re-establish agreement with the $dn_{ch}/d\eta$ distribution, as shown
in fig. 19.

  * •

StringFlav:probQQ1toQQ0join, which is really a vector of four parameters,
control the spin of the junction diquark, via a suppression of spin-1 diquark
formation relative to spin-0. This suppression factor is additional to the
factor 3 enhancement of spin-1 from counting the number of states. By
controlling spin-1 diquark production, one can directly modify probabilities
of spin-3/2 baryon production from junctions. The parameter has four
components that set the suppression factor given the heaviest quark flavour
present is $u/d$, $s$, $c$, or $b$ respectively. In CR Mode 2, these
parameters were all (arbitrarily) set to be equal to the corresponding
suppression factor in diquark string breaks, StringFlav:probQQ1toQQ0join =
{0.0275, 0.0275, 0.0275, 0.0275}, whereas the default values in Pythia are
(likewise arbitrarily) {0.5, 0.7, 0.9, 1.0}. The heavy suppression in the CR
Mode 2 values was not particularly well motivated and not tuned to any
particular set of data. Here, in view of new measurement results that have
shown comparably large $\Sigma^{0,+,++}_{c}/\Lambda_{c}^{+}$ ratios, we choose
to revert back to the default values instead of those provided by CR Mode 2.
The $\Sigma^{0,+,++}_{c}/D^{0}$ and
$\Lambda^{+}_{c}(\leftarrow\Sigma^{0,+,++}_{c})/\Lambda_{c}^{+}$ distributions
are shown in fig. 20. Though $\Sigma_{c}$ nor $\Lambda_{c}$ are spin-1/2
baryons, increasing StringFlav:probQQ1toQQ0join values enhances spin-3/2
baryon production, resulting in more excited $\Sigma_{c}$ and $\Lambda_{c}$
states. These spin-3/2 excited states decay favourably to $\Lambda_{c}$ rather
than $\Sigma_{c}$ Workman:2022ynf , and thus results in the reduction of the
$\Sigma^{0,+,++}_{c}/D^{0}$ and
$\Lambda^{+}_{c}(\leftarrow\Sigma^{0,+,++}_{c})/\Lambda_{c}^{+}$ ratios hence
describing the data more accurately. For this reason in the remainder of this
paper, we use the default Pythia values of StringFlav:probQQ1toQQ0join = {0.5,
0.7, 0.9, 1.0}. We note that these choices are obviously still to be regarded
as preliminary and subject to further smaller revisions in a more careful
tuning study.

Figure 20: $\Sigma^{0,+,++}_{c}/D^{0}$ (left panel) and
$\Lambda^{+}_{c}(\leftarrow\Sigma^{0,+,++}_{c})/\Lambda_{c}^{+}$ (right panel)
distributions from ALICE ALICE:2021rzj  data, showing the Monash (2013) tune,
and both the old (blue) and new (magenta) junction modelling. The new values
of StringFlav:probQQ1toQQ0join = {0.5, 0.7, 0.9, 1.0} are shown with solid
lines. The original CR Mode 2 values of {0.0275, 0.0275, 0.0275, 0.0275} are
shown with dashed lines.

Figure 21: The $\Lambda/K^{0}_{S}$ (top row of plots) and $\Xi/\Lambda$
(bottom row of plots) ratios. The left and right panels show these ratios as
functions of $p_{\perp}$ and $|y|$ respectively. All events are 7 TeV NSD
events from CMS CMS:2011jlm , where Pythia simulations have a lifetime cut of
$\tau_{max}=10$ mm/c and no $p_{\perp}$ cuts on final state particles.

We shall consider the impacts of the revised junction modelling on baryon-to-
meson and baryon-to-baryon ratios, first in the light-flavour sector and then
in the heavy-flavour sectors. In the light-flavour sector, we do not expect
noticeable changes, as the revised modelling should only affect a small corner
of the production phase space. Fig. 21 depicts the ratio of light-flavour
distributions, with the $\Lambda/K$ and $\Xi/\Lambda$ ratios with respect to
transverse momentum and rapidity. As expected, the new modelling is very
similar to the old, with both predicting roughly the same rapidity and
$p_{\perp}$ distributions. The $\Lambda/K$ ratio remains somewhat
overpredicted at high $p_{\perp}$, but since this is already the case for the
Monash tune it presumably originates from diquark-antidiquark string breaks
and not from junctions. Another noticeable feature of fig. 21 is the
underprediction of the $\Xi/\Lambda$ ratio, which is present in all models
shown. The $\Xi/\Lambda$ ratio is a baryon-to-baryon ratio of a double-strange
to single-strange baryon, thus this underprediction appears indicative of a
need for strangeness enhancement which cannot be described by the inclusion of
junctions alone. We plan to return to the question of strangeness enhancement
in a separate study.

Figure 22: The top row shows the $p_{\perp}$ distributions of baryon-to-meson
ratios, with the left panel showing the prompt $\Lambda_{c}^{+}/D^{0}$ ratio
from the ALICE collaboration ALICE:2021rzj  and the right panel shows the
$\Lambda_{b}^{0}/B^{0}$ from the LHCb collaboration LHCb:2023wbo . Both sets
of data are for $\sqrt{s}=13$ TeV inelastic events, with rapidity ranges
$|y|<0.5$ and $2<|y|<4.5$ respectively. The bottom row of plots shows the
$\Lambda_{b}$ asymmetry LHCb:2021xyh  for $\sqrt{s}=7$ TeV events as a
function of $p_{\perp}$ (left panel) and $y$ (right panel) in the rapidity
range $2.15<y<4.10$ and transverse momentum range $2<pT<27$ GeV.

Turning now to heavy-flavour baryon-to-meson ratios, fig. 22 shows the prompt
$\Lambda_{c}^{+}/D^{0}$ and $\Lambda_{b}^{0}/B^{0}$ ratios as a function of
$p_{\perp}$. When examining the success of the CR models, the
$\Lambda_{c}^{+}/D^{0}$ ratio has been shown to fit the data fairly well,
however in much previous literature the quality of the description of the
$\Lambda_{b}^{0}/B^{0}$ ratio had not been studied. Given the new junction
modelling, the $\Lambda_{c}^{+}/D^{0}$ is largely unchanged and exhibits
similar behaviour to the old model, with perhaps a slightly higher relative
enhancement at very low $p_{\perp}$. (We note that the old-model distributions
here are slightly lower than those typically shown in other literature such as
ALICE:2021rzj  due to the modified values of StringFlav:probQQ1toQQ0join.)
However the $\Lambda_{b}^{0}/B^{0}$ ratio is overpredicted by a factor of two
in both the old and new junction modelling. We regard this as the main issue
that remains to be addressed and would like to follow up with a dedicated
study of CR (and related) effects in $b$ vs. $c$ production. At this point, we
restrict ourselves to the following remarks: though these ratios are for
different rapidity ranges, the distributions according to Pythia are largely
unaffected by the the chosen rapidity range. Interestingly the data shows both
ratios follow roughly the same $p_{\perp}$ dependence, starting at values of
around 0.5 at low $p_{\perp}$, which neither model appears to be able to
replicate. The new junction model goes in the right direction (but not far
enough), resulting in a slightly lower $\Lambda_{b}^{0}/B^{0}$ ratio compared
to the old junction model. This is not surprising as per the results in fig.
16, which show the new modelling results in slightly harder b-baryons from
junctions, which in turn would result in slightly less soft $\Lambda_{b}$
baryons. From these distributions alone however it is unclear whether this
discrepancy between the model and data is from the overprediction of
$\Lambda_{b}$ production or underprediction of $B^{0}$ production. It would be
particularly insightful to compare to experimental results of other ratios for
these heavy-flavour hadrons, such as a $\Lambda_{b}/\Lambda_{c}$ ratio or
$B^{0}/D^{0}$ ratio, to help identify where the model predictions are falling
short.

We also validate that the new modelling well predicts the overall behaviour of
the $\Lambda_{b}$ asymmetry LHCb:2021xyh , which the QCD-CR model in Pythia
has already been shown to predict quite well. This is evident in the bottom
row of plots in fig. 22 with respect to $p_{\perp}$ and $y$, and shows the
dilution of the $\Lambda_{b}$ asymmetry particularly at low $p_{\perp}$ from
junction baryons. If anything, the dilution effect seems to be somewhat
overpredicted, consistent with what we see in the $\Lambda_{b}^{0}/B^{0}$
distribution.

Figure 23: $\Xi_{c}$-to-hadron ratios as a function of transverse momentum for
$\sqrt{s}=7$ TeV events at ALICE ALICE:2021bli  at midrapidity $|y|<0.5$. Top
left: $\Xi_{c}^{0}/D^{0}$. Top right: $\Xi_{c}^{+}/D^{0}$. Bottom left:
$\Xi_{c}^{0}/\Lambda_{c}^{+}$. Bottom right:
$\Xi_{c}^{0,+}/\Sigma_{c}^{0,+,++}$.

Though junction formation provides enhancement of low-$p_{\perp}$ baryon
production, as mentioned above in relation to the $\Xi/\Lambda$ ratio fig. 21,
it is insufficient to describe strange baryon production. This is further
evident given the distributions shown in fig. 23, which depicts ratios
involving the single strange $\Xi_{c}$ baryon to non-strange hadrons. One can
clearly observe that Pythia underpredicts these sets of data by a significant
amount, regardless of which form of junction modelling is used, demonstrating
that additional baryon enhancement from junctions is insufficient to describe
strange baryon production rates. Alongside other strange-hadron data such as
the strange hadron-to-pion ratios from ALICE ALICE:2016fzo , there is a clear
indication for a need for an additional strangeness enhancement mechanism. An
example of such a strangeness-enhancement model which is already implemented
in Pythia would be the Rope hadronization mode Bierlich:2016faw ;
Bierlich:2017sxk . Though we will not elaborate on strangeness further in this
paper, we include these distributions here to emphasise that CR and the
enhanced baryon production from junctions alone are not sufficient to describe
strange-particle distributions.

##  6 Conclusions and Outlook

We have revisited the equations of motion for string junctions with special
attention to limits of soft and/or slow endpoints, for which the previous
fragmentation framework in Pythia was not robust. Based on our considerations,
we have developed a revised fragmentation of string junctions which we have
implemented in Pythia. A main part of this work has been an updated procedure
to find the so-called “average” junction rest frame. This was largely
motivated by the lack of special treatment of soft-leg junctions as well as
the previous procedure resulting in convergence errors for around 10% of
events.

In the updated modelling, the reliance on an iterative procedure has been
removed entirely. We do so by stepping away from the notion that the average
JRF must be some overall Mercedes frame (which does not exist in all cases),
and instead map out the junction motion over time and then explicitly average
this behaviour. We also ensure to include special treatment of soft junction
legs, including the mapping out of oscillations of a soft leg around a
junction and the novel pearl-on-a-string treatment. Though we do not expect
large differences in hadron distributions given the updated modelling, the
overall treatment of junctions outlined in this paper results in a more stable
and, we believe, reliable junction hadronization procedure that is founded in
more rigorous modelling. This allows us to draw stronger conclusions on
junction predictions. For future studies a full retuning effort should be
undertaken, including both these junction revisions along with all other code
updates in the updated version of Pythia 8.311.

Further experimental data for comparisons would be of great interest. Most
interesting perhaps would be further studies vs. multiplicity measures, as
well as heavy-flavour hadron ratios, such as $\Lambda_{b}/\Lambda_{c}$ and
$B^{0}/D^{0}$ ratios, and studies of global event properties in events that
contain low-$p_{\perp}$ charm or beauty hadrons. This could assist in finding
the source of the discrepancy between the string modelling and experimental
results given the current overprediction of the $\Lambda_{b}/B^{0}$ ratio.
Further studies of baryon species sensitive to the spin-1 vs. spin-0 diquark
compositions would also be revealing, and e.g. extending this to the $b$
sector; as would studies of multiply-heavy-flavoured baryons. Baryon
correlations would also be very insightful as it may give some indication to
the production source of baryons. Diquark-antidiquark pair creation via
standard string breaks results in a baryon produced next to an antibaryon,
resulting in strong correlations. The popcorn mechanism Andersson:1984af ;
Eden:1996xi  for diquark production allows for meson production between a
diquark-antidiquark pair, resulting in more mild baryon-antibaryon
correlations. Junction baryons are expected to be even less correlated.

There is also reason to reexamine $e^{+}e^{-}\rightarrow WW$ events with
junctions in view of the possible formation of junction-antijunction
configurations there. An example of such a process is depicted in the top
image of fig. 24, with the bottom image of fig. 24 showing the two possible
string configurations for such a process; the LC dipole strings and the
possible beyond-LC junction-antijunction configuration. These junction-
antijunction topologies are formed when there are two quarks moving in one
direction, and two antiquarks moving in the opposite direction. As a result,
one would expect to find an enhancement of baryons, originating from a
junction baryon and antibaryon moving in opposite directions, a behaviour not
captured by previous $e^{+}e^{-}$ models. Measurements of possible baryon
enhancement from these junction-antijunction structures along with the
correlation of such baryons would be of particular interest for future studies
into junction modelling. Additionally there are many other reasons to further
explore $e^{+}e^{-}$ collision systems as several interesting experimental
observables were not considered during the LEP experiment such as multiplicity
dependent measurements, cf., e.g., Hunt-Smith:2020lul . These could also form
targets for fragmentation studies at future $e^{+}e^{-}$ colliders, where even
relatively rare (e.g., multi-strange, heavy-strange, etc.) baryon species can
be probed with high statistics.

Particularly careful examination of the density of string systems is also
required for future studies. As mentioned above, the inclusion of junctions
alone is insufficient to describe the strange $\Xi$ and $\Xi_{c}$ hadron
ratios, which along with multiplicity-dependent light strange-hadron-to-pion
ratios from ALICE ALICE:2017jyt  is indicative of a need for a strangeness-
enhancement mechanism that scales with system size. Though existing proposals
for such mechanisms such as rope hadronization Bierlich:2017sxk  have already
shown the ability to describe strange-hadron-to-pion ratios from ALICE, they
appear to remain insufficient in describing the enhancement seen in $\Xi_{c}$.
Further, baryon enhancement at low $p_{\perp}$ does not seem the solution
given the well described $\Lambda_{c}/D^{0}$ ratio and the overpredicted
$\Lambda_{b}/B^{0}$ ratio. Hence one must question whether strangeness
enhancement for heavy-flavour baryons occurs in some subtly different way than
for the light flavours. It would also be insightful to examine $\Xi_{b}$
hadron ratios, particularly comparison of the $\Xi_{b}/\Lambda_{b}$ to
$\Xi_{c}/\Lambda_{c}$ or $\Xi_{b}/B^{0}$ to $\Xi_{c}/D^{0}$ ratios, which may
help pick apart the patterns of baryon- vs. strangeness-enhancement
mechanisms.

In fact, our Pythia 8.311 implementation also includes the ability to modify
the strangeness probability for breaks around junctions themselves, which may
in turn result in more heavy strange baryons in particular and perhaps assist
in describing the $\Xi_{c}$ baryon distribution. A detailed examination into
such modelling however is left to a future study.

More generally, further studies into the string environment dependence of
fragmentation processes would be interesting, particularly the examination of
differences between triplet and octet string fragmentation. Experimentally
this could correspond to looking at the tips of gluon jets, LEP data of the
so-called hairpin configuration, and perhaps even diffractive LHC events
(which may be enriched in gluon fragmentation without significant transverse
excitations) may also be able to provide useful insight.

Figure 24: Top panel: $e^{+}e^{-}\rightarrow WW$ process. Both $W$ bosons
decay to a $q\bar{q}$ pair with colours green-antigreen and red-antired. In
this example the two (anti)quarks move are moving in the same direction.
Bottom panel: the left image shows the LC string configuration given the
$e^{+}e^{-}\rightarrow WW$ process above, and the right image shows a possible
junction-antijunction string configuration formed after colour reconnections
given the red-green combination results in an antitriplet.

There are of course also many more general open-ended questions within the
modelling of strings which one can tackle, including alternatives to the
Schwinger mechanism (e.g., Fischer:2016zzs ) and looking beyond the Lund
string model. The Lund model assumes a constant string tension motivated by
the linear portion of the QCD Cornell potential. Lattice support for this is
based on assuming static colour charges, and the evidence from hadron
spectroscopy also has elements of steady-state evolution. However, in a
hadronization model we are concerned with dynamic strings which are expanding
at near the speed of light and fragmenting as they do so. This begs the
question whether a straightforward model with a constant string tension is
sufficient. Few forays into more general models have so far been done (e.g.,
Hunt-Smith:2020lul ) and more would certainly provide for interesting future
extensions to the string modelling in Pythia.

###### Acknowledgements.

We are grateful to Gösta Gustafson and to Torbjörn Sjöstrand for many
insightful discussions of the motion and fragmentation of junction string
systems. We also thank the staff and students at the Rudolf Peierls Centre for
Theoretical Physics in Oxford for providing a wonderful working environment
during the completion of this work. This work was funded by the Australian
Research Council via Discovery Project DP230103014 “Beautiful strings”, by the
Royal Society Wolfson Visiting Fellowship “Piercing the precision barrier in
high-energy particle physics”, and by the Monash-Warwick Alliance for Particle
Physics. The work was also supported in part by the European Union’s Horizon
2020 research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 722105 — MCnetITN3.

##  Appendix A Model Parameters

The parameters used in the models in this paper (specifically the ones that
differ from the Monash tune and CR Mode 2) are listed in tab. 1.

|  Monash tune |  “Old” model Pythia 8.310 |  “New” model Pythia 8.311 [gluon approximation ]
---|---|---|---
StringPT:sigma | = 0.335 | = 0.335 | = 0.335
StringZ:aLund | = 0.68 | = 0.36 | = 0.36
StringZ:bLund | = 0.98 | = 0.56 | = 0.56
StringZ:rFactC | = 1.5 | = 1.5 | = 1.5
StringFlav:mesonCvector | = 1.5 | = 1.5 | = 1.5
StringFlav:probStoUD | = 0.217 | = 0.2 | = 0.2
StringFlav:probQQtoQ | = 0.081 | = 0.078 | = 0.078
StringFlav:probQQ1toQQ0join | = 0.5, | = 0.5 | = 0.5
|  0.7, |  0.7, |  0.7,
|  0.9, |  0.9, |  0.9,
|  1.0 |  1.0 |  1.0
MultiPartonInteractions:pT0Ref | = 2.28 | = 2.25 | = 2.25
BeamRemnants:remnantMode | = 0 | = 1 | = 1
BeamRemnants:saturation | – | = 5 | = 5
ColourReconnection:mode | = 0 | = 1 | = 1
ColourReconnection:allowDoubleJunRem | = on | = off | = off
ColourReconnection:m0 | – | = 0.3 | = 0.3
ColourReconnection:allowJunctions | – | = on | = on
ColourReconnection:junctionCorrection | – | = 1.20 | = 1.20
ColourReconnection:timeDilationMode | – | = 2 | = 2
ColourReconnection:timeDilationPar | – | = 0.18 | = 0.18
ColourReconnection:lambdaForm | = 0 | = 0 | = 0
StringFragmentation:pNormJunction | = 2 | = 2 | = 2
StringFragmentation:pearlFragmentation | = off | = off | = off [= on]
StringFragmentation:pearlProbFactor | – | – | – [= 4]
Table 1: Parameter values used for the model comparisons in this study.

##  Appendix B Pearl-on-a-string derivation

Here we show the derivation of eq. (15) given the Ariadne frame setup with two
back-to-back massless quarks aligned along the z-axis and a massive quark with
momentum along the x-axis, with pearl-on-a-string cases corresponding to the
massive quark has having initial momentum $p_{0}$ such that the initial
velocity is less than 1/2.

Figure 25: Single string piece near the massive pearl quark at times $t$ and
$t+dt$, with transverse velocity $v_{\perp}$, and string length $dl$ gained
due to the motion of the pearl quark.

We study the motion of the massive quark by looking at a string pieces close
to the pearl quark at times $t$ and $t+dt$ as depicted in fig. 25. By