# Simulating Glueball Production in $N_{f}=0$ QCD

David Curtin<EMAIL_ADDRESS>Department of Physics, University of
Toronto, Canada  Caleb Gemmell<EMAIL_ADDRESS>Department of
Physics, University of Toronto, Canada  Christopher B. Verhaaren
<EMAIL_ADDRESS>Department of Physics and Astronomy, Brigham Young
University, Provo, UT 84602, USA

###### Abstract

In an $SU(N_{c})$ gauge theory with zero light quark flavours $N_{f}=0$, the
only hadronic states that form below the confinement scale are composite gluon
states called glueballs. These minimal confining sectors arise in many Hidden
Valley extensions of the Standard Model, including scenarios that could hold
the solution to the dark matter question and the hierarchy problem.
Quantitative study of dark glueball phenomenology requires an understanding of
pure glue hadronization, which to date is severely lacking. In this work we
show that significant progress can be made by combining a perturbative pure
glue parton shower with a self-consistent and physically motivated
parameterization of the unknown non-perturbative physics, thanks to the modest
hierarchy between the glueball mass and the confinement scale. We make our
simulation code available as the public GlueShower package, the first glueball
generator for Hidden Valley theories, and perform preliminary studies of
several glueball production observables, with theoretical uncertainties that
take the full range of possible hadronization scenarios into account. We hope
this will enable new studies of dark sector phenomenology that were previously
inaccessible.

###### pacs:

Valid PACS appear here

††preprint: PREPRINT

##  I Introduction

The ongoing mysteries of the nature of dark matter (DM) and the electroweak
hierarchy problem have long been driving forces for extending the Standard
Model (SM). The hierarchy problem motivates searches for new particles below
the TeV scale but current experiments have not yet found evidence for their
existence. To solve these problems, increasingly interesting and complex dark
sectors are being considered. A popular framework that aims to address these
issues are Hidden Valley (HV) models Strassler and Zurek (2007). These often
include SM-singlet particles charged under a confining $SU(N_{c})$ group, see
e.g. Kang and Luty (2009); Bai and Schwaller (2014); Renner and Schwaller
(2018); Mies  _et al._ (2021). Couplings to the SM are possible via portal
interactions Holdom (1986); Patt and Wilczek (2006); Falkowski  _et al._
(2009) which tend to be very weak, allowing the possibility of GeV scale
states that can evade current experimental bounds.

Some realizations of the HV framework, generally referred to as neutral
naturalness models, solve the little hierarchy problem by cancelling quadratic
SM contributions to the Higgs mass with particles uncharged under SM color.
Important examples include Mirror Twin Higgs Chacko  _et al._ (2006a),
Fraternal Twin Higgs Craig  _et al._ (2015), Folded Supersymmetry Burdman  _et
al._ (2007), and many more Barbieri  _et al._ (2005); Chacko  _et al._
(2006b); Cai  _et al._ (2009); Poland and Thaler (2008); Cohen  _et al._
(2018); Cheng  _et al._ (2018). Additionally, HV models can produce rich and
diverse phenomenology, distinct from any SM processes, such as long lived
particles (LLPs) Alimena  _et al._ (2020); Aad  _et al._ (2013, 2020), soft
unclustered energy patterns (SUEPs) Strassler (2008); Knapen  _et al._ (2017);
Barron  _et al._ (2021), and ‘dark showers’ Knapen  _et al._ (2021), leading
to semi-visble or emerging jet signatures Cohen  _et al._ (2015, 2017, 2020);
Schwaller  _et al._ (2015); Linthorne and Stolarski (2021); Sirunyan  _et al._
(2019); Tumasyan  _et al._ (2021).

An important special case of Hidden Valleys is the pure Yang-Mills QCD case
with $N_{f}=0$. The only hadronic states in the dark sector are a spectrum of
dark glueballs Morningstar and Peardon (1999); Lucini and Moraitis (2008);
Teper (1998); Lucini  _et al._ (2010); Athenodorou and Teper (2021); Yamanaka
_et al._ (2021a), which can decay to SM states via dimension 6 or 8 operators
Juknevich (2010); Juknevich  _et al._ (2009) and have potentially long
lifetimes on collider or even cosmological scales. $N_{f}=0$ QCD-like sectors
appear commonly in neutral naturalness models; for example, in the Fraternal
Twin Higgs Craig  _et al._ (2015) only the third generation of SM fermions is
mirrored in the dark sector, leaving no strongly interacting states below the
confinement scale. Dark glueballs can then be the lightest hadronic states in
the twin spectrum. Dark glueballs have also been considered as potential DM
candidates Faraggi and Pospelov (2002); Boddy  _et al._ (2014a, b); Garcia
Garcia  _et al._ (2015); Soni and Zhang (2016); Soni  _et al._ (2017);
Forestell  _et al._ (2017, 2018); Yamanaka  _et al._ (2021b, 2020); Jo  _et
al._ (2021), with their relatively strong self-interaction giving rise to
interesting astrophysical signatures Spergel and Steinhardt (2000); Weinberg
_et al._ (2015).

Clearly, studying this scenario in detail is highly motivated. However, to
date there is no reliably way of simulating dark glueball production in high-
energy processes. Previous studies have resorted to making very simplistic
conservative assumptions, like assuming exotic Higgs decays to just two mirror
glueballs in studies of LLP signals in neutral naturalness Curtin and
Verhaaren (2015); Chacko  _et al._ (2016). Another approach is the use of
analytical approximations for the final glueball distributions after dark
hadronization Burdman and Lichtenstein (2018), but this involves some ad-hoc
parameter choices as well as being inconsistent for all but very high initial
energies, due to the relatively high mass of glueballs compared to the
confinement scale. Clearly, the absence of a reliable event generator for
$N_{f}=0$ Hidden Valleys severely hampers their phenomenological and
experimental study Knapen  _et al._ (2021). Our work addresses this
shortcoming, opening the door to a large variety of new and detailed
investigations.

The difficulty in simulating $N_{f}=0$ QCD arises from the unknown nature of
hadronization without light quarks. For $N_{f}>0$ with some dark quark masses
below the dark confinement scale, the Lund String model Andersson  _et al._
(1983) can in principle be used to describe hadronization in the dark sector,
with the existence of light colored states allowing tubes of color flux to
break via light quark pair production.111Alternative hadronization schemes
such as the cluster model Webber (1984) implemented by HERWIG++ Bahr  _et al._
(2008) also only apply for the case with light quarks, as does preconfinement
Amati and Veneziano (1979). This is implemented in the Hidden Valley module
Carloni and Sjostrand (2010); Carloni  _et al._ (2011) of the PYTHIA 8
Bierlich  _et al._ (2022) generator. In pure $SU(N_{c})$ Yang-Mills theory, on
the other hand, no existing hadronization model has been implemented so far.

In this work we present a simulation strategy for obtaining dark glueball
final states from pair produced dark gluons with some initial center-of-mass
energy $M$, implemented as the public python package
GlueShower222github.com/davidrcurtin/GlueShower for $N_{f}=0$ $SU(N_{c})$ QCD
with $N_{c}=\\{2,\ldots,12\\}$. Because dark gluon production and decay are
highly dependent on the specific BSM model, we do not specify those aspects of
the dark shower process, instead focusing on the perturbative dark gluon
shower and hadronization into dark glueballs. This can be combined with other
event generators for production and decay to give a complete signal simulation
for a given dark sector.333A follow-up paper applying this work to study the
possible indirect detection signals from Dark Matter annihilating to dark
glueballs in our galaxy is currently in progress Curtin and Gemmell .

Apart from the practical usefulness of assembling a useable event generator
for $N_{f}=0$ QCD, the main novelty of GlueShower is our parameterization of
different possible hadronization mechanisms, given that the underlying non-
perturbative dynamics are even less well understood from first principles than
hadronization with light quarks. Our hadronization model is simple and
physically motivated, with enough built-in variation to span the space of
reasonably possible ‘jet-like’ and ‘plasma-like’ final outcomes. Despite
faithfully incorporating our large theoretical ignorance of pure Yang-Mills
hadronization, the resulting predictions are of sufficient precision to make
them highly useful for dark sector searches and constraints.

In Sec. II, we briefly review the known properties of glueballs as obtained
from lattice QCD. (For simplicity, when referring to gluon and glueballs in
this paper we refer to the pure Yang-Mills case, explicitly specifying when we
instead refer to SM gluons or states.) The perturbative aspects of QCD
relevant to our Monte Carlo generator are reviewed very briefly for
completeness in Sec. III. Section IV covers the hadronization process of
GlueShower. In Sec. V we simulate glueball production for a variety of
hadronization assumptions, define a set of 8 benchmark hadronization
parameters to span the range of physically reasonable possible outcomes, and
make new predictions for observables of glueball production with theoretical
uncertainties derived by the variation across these hadronization benchmarks.
We conclude in Sec. VI.

##  II Dark Glueballs

The properties of $SU(N_{c})$ glueballs have been studied on the lattice for
decades Morningstar and Peardon (1999); Lucini and Moraitis (2008); Teper
(1998); Lucini  _et al._ (2010); Athenodorou and Teper (2021); Yamanaka  _et
al._ (2021a), establishing a spectrum of twelve stable states in the absence
of external couplings, as shown in Figure 1. These states are distinguished by
their $J^{PC}$ quantum numbers, and their masses can be parameterised entirely
in terms of the confinement scale, $\Lambda$.

Figure 1: Glueball mass $m_{G}$ spectrum for pure $SU(3)$ Yang-Mills theory
Morningstar and Peardon (1999) in terms of the lightest glueball mass $m_{0}$,
plot taken from Juknevich (2010).

In this work we use the lattice values calculated in Athenodorou and Teper
(2021). Across the range of $N_{c}=\\{2,\ldots,12\\}$ values we consider, the
lightest glueball mass $m_{0}$ is approximately 6$\Lambda$. This is much
heavier than e.g. SM hadrons, and important for our characterization of
glueball hadronization. Additionally they also provide the glueball masses for
the next nine states and relate the string tension $\sqrt{\sigma}$ to the
three-loop confinement scale in the $\overline{MS}$ scheme. Since
$\Lambda_{\overline{MS}}$ is quoted to 3-loop accuracy we also use the 3-loop
running coupling as given in Prosperi  _et al._ (2007), which has been used
with tree-level splitting functions in other Monte Carlo generators Bahr  _et
al._ (2008).444Using the 3-loop running coupling is simply to be consistent
with the lattice results and not indicative of a higher level of accuracy, as
any improvement compared to the one-loop result is subdominant to the
uncertainties of our hadronization model. These values are provided for
$N_{c}$ between 2 and 12, and because the perturbative QCD calculations are
easily adjusted for a general $N_{c}$, GlueShower generates glueball showers
for this same range of $SU(N_{c})$ confining sectors. This opens the door for
simulating the phenomenology of a wide range of interesting exotic dark
sectors, and hopefully encourages study of confining dark sectors beyond the
SM-like $SU(3)$ case Boddy  _et al._ (2014a); Soni and Zhang (2016); Batell
_et al._ (2020); Kilic  _et al._ (2021).

Glueball wavefunctions have also been studied on the lattice and by other
methods. These have determined the average size of glueballs, with the
$r_{\text{rms}}$ value for the $0^{++}$ glueball typically found to be of
order $\sim\Lambda^{-1}$ Hou  _et al._ (2001); Ishii  _et al._ (2001); Loan
and Ying (2006). The next heaviest state, $2^{++}$, is approximately twice the
size of the $0^{++}$ state Loan and Ying (2006). There is still uncertainty in
these measurements and questions regarding whether quenched QCD lattice
studies will agree with the pure Yang-Mills results. Still, these small radii
support focussing on local physics rather than considering glueballs as truly
extended objects, which guides our discussion of hadronization below.

##  III Perturbative Shower Review

Perturbative QCD is an extensively studied and established field within
quantum field theory. Following Ellis  _et al._ (2011); Sjostrand  _et al._
(2006) we briefly review the salient details of perturbative parton showers in
this section for completeness and to establish notation. More complete
documentation of how this is implemented in GlueShower is found in Appendix A.

For $N_{f}=0$ HV models, the only parton is the gluon. Our code works in the
centre-of-mass frame of a two-gluon initial state with invariant mass
$M$.555This is the only initial state we consider, as it is of the most use to
BSM physics studies. We leave the study of other interesting cases, such as a
gas of gluons that is cooled until it undergoes confinement, for future
investigations. This amounts to simulating gluon production from the decay of
a massive scalar particle of mass $M$, but can be easily generalized to, for
example, direct di-gluon production via effective operators or intermediate
states. These initial gluons are produced with large virtualities (effective
mass-squared) $t$, and as they split into more gluons produce an increasing
ensemble of lower virtuality gluons. The energies and virtualities of these
evolving gluons are described by perturbative QCD and can be simulated as a
parton shower.

The probability that a gluon splits into two gluons with energies $z$ and
$(1-z)$, where $z$ is the energy fraction of the mother dark gluon, is
determined only by the gluon-to-gluon splitting function:

$P_{gg}(z)=2C_{A}\bigg{[}\frac{z}{1-z}+\frac{1-z}{z}+z(1-z)\bigg{]}~{},$
(III.1)

where $C_{A}=N_{c}$. This splitting function is also used to define the
Sudakov form factor, which gives the probability that the gluon evolves from
an initial virtuality, $t_{0}$, to a lower virtuality, $t$, without splitting:

$\Delta(t)=\text{exp}\bigg{[}-\int_{t_{0}}^{t}\frac{dt^{\prime}}{t^{\prime}}\int
dz\frac{\alpha_{s}}{2\pi}P_{gg}(z)\bigg{]}~{}.$ (III.2)

Note that this Sudakov form factor only accounts for the leading-log collinear
gluon enhancements, which is unable to reproduce the correct scaling of parton
multiplicity with centre of mass energy. However, this is easily rectified in
standard implementations of parton showers by imposing angular ordering on
subsequent splittings, which accounts for soft gluon interference effects at
leading order in $1/N_{c}^{2}$ Marchesini and Webber (1984), see Appendix A
for further details.

Monte Carlo parton shower evolution can be framed in the following way: Given
a gluon with virtuality and energy, $(t_{1},z_{1})$, after some step in the
evolution, what is its new virtuality and energy, $(t_{2},z_{2})$? Note that a
gluon can only decrease its virtuality by splitting. Thus, finding a value for
$t_{2}$ implies the gluon split in the intermediate step. To calculate $t_{2}$
a random number $\mathscr{R}\in[0,1]$ is generated and $t_{2}$ is found by
solving

$\Delta(t_{2})=\frac{\Delta(t_{1})}{\mathscr{R}}~{}.$ (III.3)

If there is no solution, then the gluon does not split. In this case the gluon
can only hadronize, and its shower terminates. If a $t_{2}$ is found, $z_{2}$
is determined by generating another random number
$\mathscr{R}^{\prime}\in[0,1]$ and solving:

$\int_{z_{\text{min}}}^{z}dz^{\prime}\frac{\alpha_{s}}{2\pi}P_{gg}(z^{\prime})=\mathscr{R}^{\prime}\int_{z_{\text{min}}}^{1-z_{\text{min}}}dz^{\prime}\frac{\alpha_{s}}{2\pi}P_{gg}(z^{\prime}),$
(III.4)

where $z_{\text{min}}$ is set by kinematic thresholds of the possible gluon
splittings. This Monte Carlo generation is implemented in GlueShower as
detailed in Appendix A. We also note that since $\alpha_{s}\propto N_{c}^{-1}$
at 1-loop order, and the perturbative shower is only dependent on
$\alpha_{s}P_{gg}(z)$, changes to $N_{c}$ only affect the shower due to the
slight dependence of glueball masses on $N_{c}$. In our numerical studies
below we therefore only show the $N_{c}=3$ case, with other numbers of colours
giving similar behaviour.

Lastly, as is standard in parton shower MC algorithms, we work in the leading
colour limit, which is equivalent to taking $N_{c}\rightarrow\infty$ with
$\alpha_{s}N_{c}$ kept constant ’t Hooft (1974). This amounts to ignoring
higher order colour interference effects in $1/N_{c}^{2}$ expansions.
Additionally in this limit we can use the simple representation of gluons
being the direct product of a fundamental and anti-fundamental, since the
weight of the singlet in $N_{c}\otimes\overline{N_{c}}=(N_{c}^{2}-1)\oplus 1$
vanishes. It is in this limit that the t’ Hooft double line notation can be
used to trace colour flow, as shown in Fig. 2. From SM QCD we know this
approximation works well for $N_{c}=3$, and thus also for $N_{c}>3$. It is
possible that this approximation introduces larger errors for $N_{c}=2$, but
we leave the problem of including the subleading colour corrections for future
investigations.

##  IV Hadronization

Figure 2: Comparison between Feynman diagram-like (left) and flux string-like
(right) cartoons of the evolution of colour singlet separation in $N_{f}=0$
QCD.

While the perturbative aspects of the parton shower have been extensively
studied, the non-perturbative aspects of QCD have not been computed from first
principles. We follow the usual approach, outlined in the previous section,
whereby we iteratively evolve gluons from their initial virtuality and undergo
splittings above some hadronization scale $\Lambda_{\text{had}}$, which is
typically taken to be of the same order as the confinement scale $\Lambda$
(though it turns out this assumption must be modified for our case). Some
hadronization model capturing the non-perturbative dynamics of glueball
formation is then needed to turn the gluon final states of the shower into
physical final states. In SM QCD, phenomenological hadronization models
Andersson  _et al._ (1983); Webber (1984) can be tuned to reproduce the
observed data, but these rely on the existence of light quarks, which either
allow a flux tube string to fragment, or allow colour singlet clusters to form
during preconfinement Amati and Veneziano (1979). For pure Yang-Mills, the
corresponding hadronization models have not been formulated and at any rate
cannot be tuned to data prior to the discovery of a corresponding dark
sector.666A highly exciting possibility is understanding pure glue
hadronization from first principles on the lattice. This is very challenging,
but there has been significant recent progress for $N_{f}>0$, see e.g. Bulava
_et al._ (2019). We therefore provide a novel parameterization of the various
possibilities for $N_{f}=0$ QCD hadronization, relying on simple physically
motivated arguments to capture the range of both ‘jet-like’ or ‘plasma-like’
possibilities for how pure glue might hadronize, producing a consistent
spectrum of simulated glueball final states obeying full energy-momentum
conservation. This accurately captures present theoretical uncertainties while
allowing for the quantitative study of pure glue hidden valleys with
sufficient precision for searches and constraints.

###  IV.1 An Intuitive Cartoon

We begin with a representative cartoon to guide our intuition for pure glue
hadronization. In this cartoon we assume there is only one glueball species
and consider the simple example of producing two gluons that eventually
hadronize into two glueballs. Because the glueballs are colour-singlet final
states, and we begin with two back-to-back colour-octet gluons, long-range
colour exchange must occur at some point along the shower/hadronization _if
the center-of-mass energy is high enough_.777The minimum initial energy
required for separation into multiple colour singlets will be made more
precise in the next subsection. Thus, before the glueballs can form, at
minimum one of the original gluons must split into three gluons, with one new
gluon joining the other branch and allowing colour singlet states to form,
referred to as colour octet neutralisation Minkowski and Ochs (2000).888Other
detailed splitting histories are of course possible but do not change our
qualitative argument. This illustrative cartoon is depicted in a Feynman-
diagram-like way on the left side of Fig. 2, showing gluon splittings and flow
of colour charge.

Alternatively, one can consider how tubes of colour flux in the fundamental
representation evolve throughout the shower. This equivalent representation is
depicted on the right side of Fig. 2. The initial state is an overall colour
singlet loop, with the gluons behaving as localised energy or ‘kinks’ in the
loop. Forming two final state glueballs can only be accomplished by the
initial loop fragmenting in two. As the loop evolves, gluon splitting
introduces new kinks and new colors for the flux tube segments in the loop.
Towards the end of the shower, step 3, enough gluons have formed to allow two
flux tube segments in the loop to have the same colour, making fragmentation
into two colour singlet glueballs possible by crossing the same-colour flux-
tube segments. _Thus, we see that glueball hadronization can be pictured as
crossing color-fundamental flux tubes ‘pinching off’ one flux tube loop into
two separate loops, which becomes possible even after a small number of gluon
splittings._

The tube crossing picture is supported by the fact that loop fragmentation is
hugely energetically preferred once it is possible, since flux tube crossing
has little energy cost compared to the energy stored in flux tubes of length
$\gg\Lambda^{-1}$. Making the crossing possible only requires the exchange of
arbitrarily soft IR gluons carrying colour information, which carries no
particularly relevant suppression and is in fact enhanced by $\alpha_{s}$
running to large values at low energies. This is also consistent with glueball
radii being of order the confinement scale, suggesting their formation is
dominated by local physics.

We can compare the glueball flux tubes with the Lund String Model used for
QCD. For the simple case of quark pair production, the shower and subsequent
hadronization process essentially follows the evolution of an open flux tube
with the quarks at either end. These flux tubes are easily broken through
light quark pair production, which increases the open flux tube multiplicity.
The final state flux tube states with quarks and antiquarks on either end are
then associated with a variety of mesons. $N_{f}>0$ QCD is predominately
determined by the fragmentation of open strings, while the final state for
$N_{f}=0$ QCD HV theories is determined by closed loops pinching and
fragmenting.

In summary, formation of individual glueballs requires crossing flux tubes,
which seems to be entirely unsuppressed once a small number of (possibly very
soft) gluons have been exchanged between different branches of the shower.
This required gluon exchange in turn becomes unsuppressed once the
virtualities in the shower approach the confinement scale. _This motivates the
intuition that glueball formation should occur around the confinement scale,
qualitatively similar to SM jet-like behaviour._

###  IV.2 GlueShower Hadronization Implementation

We now describe our pure glue hadronization algorithm in detail, first for the
default jet-like assumption as motivated above, and then for more exotic
plasma-like behavior, to cover all physically reasonable possibilities.

####  IV.2.1 Jet-like glueball hadronization

The assumption that glueball production proceeds in a qualitatively similarly
jet-like fashion as for SM QCD is well-motivated by physical arguments, as
explained in the previous section. However, the detailed question still
remains how to map some arrangement of gluons at the end of the perturbative
shower to a set of final state glueballs. One could implement a full Lund-like
string model of closed flux tubes and their crossing, and vary its parameters
to obtain a range of possibilities for glueball hadronization. This would be
quite involved, and we leave this for future investigations. However, thanks
to the significant separation between the confinement scale $\Lambda$ and the
lightest glueball mass $m_{0}$, we can already make great progress with much
simpler physical arguments.

We begin by assuming that there is just a single glueball species of mass
$m_{0}$, and consider a gluon in the perturbative shower that has virtuality
$\sqrt{t}>2m_{0}$ and does not undergo further splittings above scale
$2m_{0}$. If the gluon were allowed to evolve further down in virtuality, any
subsequent splittings cannot result in more than one on-shell glueball unless
there is significant momentum exchange with other branches of the shower.
Instead, subsequent splittings result in lower-energy gluons physically
clustered around the original gluon 4-momentum in a decay cone that is at most
as wide as one arising from the decay of a particle with mass below $2m_{0}$.
As $\sqrt{t}\to\Lambda$, hadronization must combine these soft gluons (plus
some soft IR gluons to exchange colour information with other branches of the
shower) into an on-shell, colour-singlet glueball with mass $m_{0}$.

Therefore, we argue that for a given shower history, an _upper bound_ on the
number of glueball states produced can be obtained by setting the
hadronization scale to $\Lambda_{\text{had}}=2m_{0}$ and simply turning the
gluons at the end of the shower into on-shell glueballs of mass $m_{0}$. Note
that $2m_{0}\sim 12\Lambda$ is well within the perturbative regime, so the
simulated momenta of gluons at the end of the shower are highly reliable.
However, in converting these gluons to glueballs some soft gluons must be
exchanged with other branches of the shower to form colour singlets. This
suggests momentum transfers of order $\Lambda/2m_{0}$, but the relatively high
mass of glueballs makes this correction factor smaller than $\approx 10\%$.
Therefore, naively turning gluons with $\sqrt{t}=\Lambda_{\text{had}}$ into
on-shell glueballs is likely to be a good approximation.999This method brushes
over the specifics of how colour information is exchanged to create colour
singlets, reminiscent of early independent fragmentation models Field and
Feynman (1978). However, these early models proved to be historically useful
in SM QCD, are even more useful for pioneering BSM studies in our case due the
larger mass of the glueball hadrons in $N_{f}=0$ compared to the typical
momentum transfer involved in colour exchange.

How could this upper bound on glueball multiplicity be violated? First, it may
be possible for two neighbouring branches of the shower to exchange gluons
with momenta between $2m_{0}$ and $\Lambda$, allowing two gluons with
virtuality $2m_{0}$ each to result in three instead of two final state
glueballs. This can still be regarded as a (marginally) perturbative process
in the shower’s regime of validity. On the other hand, this assumption can be
violated entirely in the non-perturbative regime if the flux strings between
two branches contain enough energy to produce a third glueball.

We first discuss the perturbative possibility, where the size of the momentum
exchange gives us some hope of using the perturbative shower and simple phase
space arguments to estimate this rate $P_{2\to 3}$ of turning two gluons at
the end of the shower into three glueballs. We overestimate $P_{2\to 3}$ to
show that it is small enough to ignore at our current level of precision.

Consider two gluons 1 and 2 at the end of the shower, both with virtuality
$2m_{0}$, energies $E_{1,2}$ and angle $\theta$ between their momenta. To
overestimate $P_{2\to 3}=P_{2\to 3}(E_{1},E_{2},\theta)$, we assume that a
third glueball is formed if the gluons split $1\to 1^{\prime}3$ and $2\to
2^{\prime}4$ such that $m_{1^{\prime},2^{\prime}}>m_{0}$, the daughters have
sufficient invariant mass to form the glueball $m_{34}\geq m_{0}$ and the two
daughter momenta are close in phase space compared to the confinement scale,
i.e. $\Delta p_{34}\equiv|\vec{p}_{3}-\vec{p}_{4}|<a\Lambda$ for
$a\sim\mathcal{O}(1)$, to enable the merger. We choose $a=2$ but the precise
value does not significantly affect our result. We can therefore estimate

$\displaystyle P_{2\to 3}(E_{1},E_{2},\theta)$ $\displaystyle=$
$\displaystyle\int dt_{1}dt_{2}dz_{1}dz_{2}$ $\displaystyle P_{1\to
1^{\prime}3,2\to
2^{\prime}4,m_{1^{\prime},2^{\prime}}>m_{0}}(E_{1},E_{2},\theta)$
$\displaystyle P_{m_{34}>m_{0},\Delta
p_{34}<a\Lambda}(|p_{3}|,|p_{4}|,\theta)$

This integrates over all possible splittings $(t_{1},z_{1}),(t_{2},z_{2})$ of
the two parent gluons that produce daughters 3 and 4. The first term
(splitting probability) just evaluates the Sudakov to give the splitting
probability of both parent glueballs producing two daughters with momentum
$p_{3},p_{4}$ such that the parents still have sufficient virtuality to form
their own glueballs. The second (merger probability) term is the probability,
given a random emission angle for each daughter in the transverse plane of the
parent momentum, that the two daughters could in principle combine to form a
third glueball according to our above criteria.

To simplify evaluation of this integral, we overestimate both terms
separately. The splitting probability (first term) is overestimated by letting
both gluons 1 and 2 run from $\sqrt{t}=2m_{0}$ down to $\sqrt{t}=2\Lambda$,
allowing for emission of daughters with virtualities as low as $\Lambda$. This
is a huge overestimate since we are not enforcing
$m_{1^{\prime},2^{\prime}}>m_{0}$, and turns $P_{1\to 1^{\prime}3,2\to
2^{\prime}4,m_{1^{\prime},2^{\prime}}>m_{0}}(E_{1},E_{2},\theta;p_{3},p_{4})$
into a constant $P_{\text{split}}^{\text{max}}$ that sits outside the
integral. We then set the virtualities of the daughters $3,4$ to be the
largest possible masses allowed by a given splitting to define their 4-momenta
$p_{3},p_{4}$, which in turn maximizes the merger probability (second term).
This allows us to define

$\displaystyle P_{2\to 3}^{\text{max}}(E_{1},E_{2},\theta)$
$\displaystyle\equiv$ $\displaystyle P_{\text{split}}(E_{1},E_{2})\times$
$\displaystyle\int
dt_{1}dt_{2}dz_{1}dz_{2}P_{\text{merge}}(|p_{3}|,|p_{4}|,\theta)$

satisfying $P_{2\to 3}^{\text{max}}>P_{2\to 3}$. Having thus obtained a
function $P_{2\to 3}^{\text{max}}(E_{1},E_{2},\theta)$, it is then
straightforward to generate events using the perturbative pure-gluon shower
for a variety of initial center-of-mass energies $M$, terminate at
hadronization scale $\Lambda_{\text{had}}=2m_{0}$, and evaluate the chance of
obtaining an extra glueball compared to our naive expectation of turning
gluons with virtuality $2m_{0}$ into glueballs:

$P_{\mathrm{extra\ glueball}}=1-\prod_{i,j}(1-P_{2\to
3}^{\text{max}}(E_{i},E_{j},\theta_{ij}))$ (IV.7)

where the product is over all gluon pairs $(i,j)$ in the event. Averaged over
all events, this probability is shown in Fig. 3 for $N_{c}=3$, but the result
is nearly identical for other numbers of colours. As expected, the chance of
producing additional glueballs beyond our upper bound increases with center-
of-mass energy since this gives more gluons and more chances for the required
splittings and mergers. However, for $M\lesssim 100m_{0}$, the error
introduced by ignoring $2\to 3$ production of glueballs is smaller than 5%,
and we are justified in ignoring it for our current implementation.

We now consider non-perturbative effects that could produce additional
glueballs, orienting ourselves in the SM analogues of Schwinger pair
production in QED Schwinger (1951), or pion production in the “snapping” of
colour-fundamental strings Casher  _et al._ (1979). The latter case is of most
physical interest. The energy density of the colour string is given by the
string tension $\sigma$, and the production rate of hadrons in string
fragmentation scales as

$P_{\mathrm{hadron}}\propto\text{exp}\bigg{[}-\frac{\pi
m_{\mathrm{hadron}}^{2}}{\sigma}\bigg{]},$ (IV.8)

since the constituent masses of the partons have to be provided by the colour
background field. Ignoring any additional suppressions that may arise from the
required geometrical arrangement of two colour flux tubes to allow for the
creation of a glueball, the fact that the glueball mass is much larger than
the mass of SM pions compared to the string tension should result in a large
suppression on the number of glueballs produced from the dynamics of the non-
perturbative colour strings. For $N_{f}=0,N_{c}=3$, the ratio of the
exponential factors is

$\frac{P_{\mathrm{GB}}}{P_{\mathrm{\pi}}}\sim 10^{-16}\ .$ (IV.9)

where the numerator was evaluated using
$\Lambda_{\overline{MS}}/\sqrt{\sigma}=0.5424$ and the lightest glueball mass,
$m_{GB}=6.28\Lambda_{\overline{MS}}$ Athenodorou and Teper (2021). Note that
for different $SU(N)$ groups, the constants change slightly but the extreme
suppression persists. We therefore conclude that both perturbative and non-
perturbative effects should not invalidate our jet-like approach.

Figure 3: Upper bound on the probability of producing more glueballs of mass
$m_{0}$ in gluon pair production events than the naive upper multiplicity
bound obtained by terminating the pure glue shower at
$\Lambda_{\text{had}}=2m_{0}$, as a function of the center-of-mass energy $M$.

If setting $\Lambda_{\text{had}}=2m_{0}$ and simply turning final state gluons
into on-shell glueballs gives the largest possible glueball multiplicity per
event, how can we take into account the possibility that the real number of
produced glueballs could be lower? Physically, this would correspond to colour
singlets forming via gluon exchange at a higher virtuality scale than naively
expected, which is certainly a possibility given the unknown details of non-
perturbative $N_{f}=0$ QCD. We can obtain consistent events representing this
scenario by simply terminating the shower at a higher scale before turning
final state gluons into on-shell glueballs. In other words, we set

$\Lambda_{\text{had}}=c\times 2m_{0}$ (IV.10)

where $c>1$ is a dimensionless parameter encoding our assumption of the higher
scale $\Lambda_{\text{had}}$ where the shower fragments into disconnected
colour singlets that each yield one glueball. Varying $c\sim\mathcal{O}(1)$
gives us a controlled way to parameterize different possible assumptions on
jet-like glueball hadronization, and explore the effect of this uncertainty on
observables while still producing fully consistent shower histories that
conserve energy and momentum for arbitrary choices of $c\geq 1$ and
$M/m_{0}>2$.

So far, we have only considered a single glueball species, but in reality
there is a spectrum of roughly a dozen different glueball states with
different $J^{PC}$ quantum numbers. While our arguments support the notion
that the _inclusive_ glueball _multiplicity_ and _momentum distributions_ are
dominated by local jet-like physics, the same is not true for _exclusive_
distributions for each glueball species. Soft gluon exchange at the scale
$\sim\Lambda$ (also the scale of mass differences between glueball species)
can easily exchange angular momentum and other quantum numbers to turn a given
candidate gluon into a variety of different glueball species, and rigorously
analyzing these non-local effects is far beyond our scope. We therefore adopt
a much simpler approach of assuming that the relative multiplicities of
different glueballs follows a thermal distribution Falkowski  _et al._ (2009)
in the absence of other threshold effects, with the probability for producing
glueball state $J$ given by

$P_{J}\propto(2J+1)\bigg{(}\frac{m_{J}}{m_{0}}\bigg{)}^{3/2}e^{-(m_{J}-m_{0})/T_{\text{had}}}~{},$
(IV.11)

where the glueball masses for different $N_{c}$ are known from the lattice,
and we define a hadronization temperature

$T_{\text{had}}=d\ T_{c}$ (IV.12)

which is related by a dimensionless coefficient $d\sim\mathcal{O}(1)$ the
critical temperature of the $N_{f}=0$ QCD phase transition $T_{c}$. This is
justified by investigations of final state distribution from closed string
emission Manes (2002) that support a thermal model, where $T_{c}$ is taken to
be the Hagedorn temperature. Consequently, we assume the Hagedorn temperature
is the critical temperature of deconfinement Blanchard  _et al._ (2004);
Noronha-Hostler  _et al._ (2010). In SM QCD, the critical temperature is 150
MeV Petreczky (2012), smaller than the confinement scale. In $N_{f}=0$ Yang-
Mills theories the relation between the critical temperature and string
tension has been studied on the lattice Lucini  _et al._ (2004); Boyd  _et
al._ (1996); Lucini  _et al._ (2005, 2012), with Lucini  _et al._ (2012)
finding the relation

$\frac{T_{c}}{\sqrt{\sigma}}=0.5949+\frac{0.458}{N_{c}^{2}}~{},$ (IV.13)

giving a critical temperature that is slightly larger than the confinement
scale. By combining this result with the three-loop relation between the
confinement scale $\Lambda$ and the string tension, the relative glueball
multiplicities are entirely determined by the number of colours $N_{c}$ and
the nuisance parameter $d$.

We incorporate these probabilities into our hadronization routine in the
following manner. Once it is determined that a gluon can no longer split
during the perturbative shower, it remains with some virtuality
$\sqrt{t}\geq\Lambda_{\text{had}}=2cm_{0}$. This gluon then selects a random
glueball final state with on-shell mass below its current virtuality, weighted
by the probabilities in Eq. (IV.11). This treats the thermal probabilities
$P_{J}$ as fundamental, and introduces some additional threshold effects that
favour light glueball production, since a high-virtuality-gluon that does not
split before termination of the shower can have more glueball final states
kinematically accessible to it than a gluon that was produced during the
shower with virtuality close to $2\,c\,m_{0}$. The relative multiplicity of
glueball species produced by this hadronzation routine will therefore skew
towards lighter flavours than the thermal $\\{P_{J}\\}$ alone, but we believe
this is a physically reasonable prediction of the perturbative shower based on
little more than phase space arguments.

Additionally, colour rope hadronization Biro  _et al._ (1984) in SM QCD can
affect the final state relative multiplicities by enhancing strangeness
production Bierlich (2018) via an increased string tension. This could be
relevant for the $N_{f}=0$ QCD case since the two fundamental flux strings
between gluons would behave collectively when this effect is included,
changing the relative multiplicity of glueballs. Including these dynamics is
beyond our scope, but we account for their possible effect by not treating
$T_{\text{had}}$ as firmly determined, and vary $d$ to parameterize the theory
uncertainty in relative glueball species multiplicity. Higher hadronization
temperatures favour production of heavier glueball states. Depending on the
full details of the theory and the operator by which glueballs decay to SM
states Juknevich (2010); Juknevich  _et al._ (2009), this can significantly
affect the visible phenomenology, and it is important to treat
$d\sim\mathcal{O}(1)$ as a nuisance parameter in quantitative analyses. We
demonstrate this by studying some benchmarks below.

While we expect our approach of randomly assigning glueball species identity
based on thermal probabilities to be reasonable on average over many simulated
events, this simple approximation is unlikely to give accurate intra-event
correlations between separations in momentum space and relative species
assignments of different glueballs, for example. Even so, this simple
parameterization should be sufficient for many first phenomenological
investigations.

####  IV.2.2 Plasma-like glueball hadronization

While the jet-like showering and hadronization behaviour is highly physically
motivated, we want GlueShower to cover the largest range of physically
possible showering behavior. We therefore also consider a much more exotic
plasma-like or SUEP-like regime.

Within the jet-like assumption, the fragmentation of closed flux tubes is
viewed as an IR process that results in the immediate formation of on-shell
glueballs. However, if fragmentation into colour singlets occurs at a
significantly higher scale, then a population of high mass closed flux tubes
could be produced.101010We thank Matthew Strassler for bringing this
possibility to our attention. These states could be treated as a collection
very excited glueballs, or a hot ball of gluon plasma. Much like a quark-gluon
plasma that evaporates via quasi-isotropic emission of pions, this pure glue
plasma would evaporate by emitting glueballs approximately isotropically with
thermal momenta in its restframe.

To implement this possibility within GlueShower, we introduce boolean
parameter plasma_mode which is False by default (jet behaviour) but can be set
to True to enable plasma behaviour. In plasma mode, the shower still
terminates at a scale determined by setting $c>1$ just as in jet mode,

$\Lambda_{\text{had}}=c\times 2m_{0},$ (IV.14)

but this scale is now interpreted as the scale below which the shower
separates into singlet plasma balls, and each final-state gluon is turned into
a singlet of mass $m_{\text{plasma}}=c\,m_{0}$ instead of a glueball $J$ of
mass $m_{J}$.

The evaporation of these plasma balls into glueballs is treated analogously to
dark hadron production in SUEP scenarios Strassler (2008); Knapen  _et al._
(2017); Swisdak (2013), assuming isotropic thermal glueball emission. We
assume that the thermal glueball energy distribution and their relative
species probabilities are dictated by the same temperature
$T_{\text{had}}=d\cdot T_{c}$, see Eq. (IV.11). We borrow the SUEP-simulation
methods used in Knapen  _et al._ (2017); Swisdak (2013) to generate isotropic
glueball momenta in the restframe of each plasma. In short, isotropic glueball
momenta are successively generated with a thermal energy distribution until
additional gluon emission would cause the total invariant mass of all gluon
momenta to exceed the plasma mass $m_{\text{plasma}}$. For each glueball, its
species is randomly picked weighted by the thermal probabilities $P_{J}$ in
Eq. (IV.11). Once glueball emission is completed for a given plasma-ball, the
entire system of daughter-glueballs is slightly boosted and its kinetic
energies rescaled to exactly equalize their rest frame and invariant mass with
the original plasma ball.

In plasma mode, glueballs are therefore produced in a manner that lies between
the jet-like behavior of QCD and the purely isotropic behavior of pseudo-
conformal theories that realize the SUEP scenario: the shower generates
branches of total invariant mass $m_{\text{plasma}}=\Lambda_{\text{had}}/2$,
which become plasma balls that evaporate via isotropic thermal glueball
emission in their respective restframes.

We emphasize that this possibility is highly exotic: we regard jet-like
hadronization as far more physically motivated, since realizing plasma-like
behavior requires very unusual long-distance non-perturbative effects that
arrest further fragmentation well within the naive regime of reliability for
the perturbative shower. However, we include it in our code to make sure that
even highly exotic hypotheses for the non-perturbative dynamics of $N_{f}=0$
QCD can be qualitatively accommodated.

####  IV.2.3 Input Parameters for GlueShower

In summary, GlueShower has two purely physical parameters:
$N_{c}\in(2,3,\ldots 12)$ to specify the dark gauge group and $m_{0}$ to
specify the mass of the lightest glueball.111111The large $N_{c}$ regime can
be well-covered by simulating $N_{c}=12$. This determines the confinement
scale $\Lambda$, critical temperature $T_{c}$, and the 3-loop running coupling
$\alpha_{S}(\mu$). For a given run, one also specifies the initial center-of-
mass energy $M$ for the di-gluon initial state (corresponding to di-gluon pair
production in the decay of a scalar particle of mass $M$) as well as the
number of shower histories to simulate. Theoretical uncertainties of dark
hadronization are captured in two-and-a-half nuisance parameters: the boolean
parameter plasma_mode which specifies whether gluons hardronize in the
physically motivated jet-like or the more exotic plasma-like fashion, and
multiplicative parameters $c=\Lambda_{\text{had}}/(2m_{0})$ and
$d=T_{\text{had}}/T_{c}$ which set the hadronization/singlet formation scale
and the hadronization (and plasma-ball, if in plasma mode) temperature
respectively.

Figure 4:  Average glueball multiplicity $\langle N\rangle$ for $N_{c}=3$ as a
function of center-of-mass energy $M$ relative to the lightest glueball mass
$m_{0}$. The jet-like case with the highest possible glueball multiplicity is
shown in blue, a representative plasma-like case in orange (see legend). The
bands shows the $16\%$ to $84\%$ percentile range of multiplicities in the
event distributions for each $M/m_{0}$ (jagged since multiplicities are small
integers). The analytical expectation of Eq. (V.15), normalized to the
$M/m_{0}=100$ jet-like simulation, is shown as the dashed red line.

##  V Simulation of dark glueball final states

In this section we study events generated by GlueShower. We first explore the
basic jet-like shower case, using parameters plasma_mode = False,
$\Lambda_{\text{had}}/(2m_{0})=c=1$, and $T_{\text{had}}/T_{c}=d=1$ to
demonstrate the code and its output, as well as comparing the multiplicity
scaling and shape of fragmentation functions to analytical approximations in
their expected regime of validity. We then comment on the qualitative
differences in the plasma-like case, using parameters plasma_mode = True,
$\Lambda_{\text{had}}/(2m_{0})=4$, and $T_{\text{had}}/T_{c}=1$. Finally, we
argue that a set of 8 benchmark values for the nuisance parameters covers the
physically motivated range of hadronization possibilities for glueball
production (plasma_mode = FALSE {TRUE}, $c=1,2$ $\\{4,6\\}$, $d=1,2$) and
study the resulting range of physical predictions for some observables.

###  V.1 Jet-like Hadronization

We first demonstrate how the inclusive glueball multiplicity scales with
initial centre-of-mass energy $M$ in Fig. 4 (blue). This case of
$c=\Lambda_{\text{had}}/(2m_{0})=1$ represents the largest multiplicity
possible in the jet-like case, but even for large $M\sim 100m_{0}$, the
sizeable hierarchy between $m_{0}$ and $\Lambda_{\text{had}}$ results in only
a handfull of produced glueballs per event. This differs greatly from the high
multiplicity production of pions in high-energy QCD jets. The blue band shows
the range of multiplicities produced in simulated events, which is much more
sharply peaked than a Poisson distribution. We show the $N_{c}=3$ case but the
results are similar for other numbers of colours.

As a consistency check, we compare the average multiplicity predicted by
GlueShower to the analytical expectation for average hadron multiplicity
scaling in perturbative QCD. The standard result, for example found in Ellis
_et al._ (2011) for $N_{f}=0$, is

$\langle
N(E_{CM}^{2})\rangle\propto\text{exp}\Bigg{[}\frac{12\pi}{11C_{A}}\sqrt{\frac{2C_{A}}{\pi\alpha(E_{CM}^{2})}}+\frac{1}{4}\text{ln}\bigg{(}\alpha(E_{CM}^{2})\bigg{)}\Bigg{]}~{}.$
(V.15)

This is normalized to the absolute multiplicity of our simulation for
$M=100m_{0}$ and shown as the red dashed line in Fig. 4. Note the good
agreement for $M\gg m_{0}$,121212Note that some deviation is expected since
our shower uses 3-loop running of the coupling but as expected, the scaling
relation breaks down for $M$ closer to the glueball mass, demonstrating that
finite-mass effects make analytical approximations of glueball distributions
unreliable even for modest initial energies.

Figure 5: Glueball Energy Spectra for an example jet-like shower, with
plasma_mode = False, $\Lambda_{\text{had}}/(2m_{0})=1$,
$T_{\text{had}}/T_{c}=1$. _Left:_ Representative $0^{++}$ spectra for a range
of $M$ values. The $M/m_{0}=100$ spectrum is fit to Eq. V.16 (red) on the
range $0.1<x<1$. _Right:_ Spectra for the lightest 10 glueball states for
initial center-of-mass energy $M/m_{0}=100$.

Figure 6: Glueball Energy Spectra for an example plasma-like shower, with
plasma_mode = True, $\Lambda_{\text{had}}/(2m_{0})=4$,
$T_{\text{had}}/T_{c}=1$. _Left:_ Representative $0^{++}$ spectra for a range
of $M$ values. The $M/m_{0}=100$ spectrum is fit to Eq. V.16 (red) on the
range $0.1<x<1$. _Right:_ Spectra for the lightest 10 glueball states for
initial center-of-mass energy $M/m_{0}=100$.

We now compare the simulated glueball energy spectra to analytical
expectations. In the SM, fragmentation functions are a priori unknown, thus a
physically motivated functional form is used as an ansatz and fit to data. We
compare our jet-like output to one of the commonly used Colangelo and Nason
function Colangelo and Nason (1992)

$D_{g}^{GB}(x)\propto x^{\alpha}(1-x)^{\beta}~{},$ (V.16)

where $x=2E_{GB}/M$ and $E_{GB}$ is the glueball energy, thus $x\in[0,1]$.
This function has been used for heavy quark fragmentation Zyla  _et al._
(2020) and is equivalent to the gluon-to-kaon fragmentation function
parameterisation in de Florian  _et al._ (2017), which has also been studied
using lattice QCD Salas-Chavira  _et al._ (2021). Using a heavy quark
fragmentation function for a pure glue shower might seem surprising, but this
simply encodes that the energies of heavy hadronic final states are dominated
by the heavy quark energies. This process is determined by the perturbative
shower since it occurs significantly above the confinement scale, much like
our treatment of jet-like glueball hadronization, making such a functional
form a reasonable ansatz. A similar approach was also used to analytically
estimate dark glueball momenta in Burdman and Lichtenstein (2018).

We compare our simulated events to an analytical fragmentation function in the
high energy regime $M/m_{0}=100$ and for $0.1<x<1$ de Florian  _et al._
(2017), where finite mass effects are less significant. As shown in Fig. 5
(left), we find very good agreement in this regime of applicability, provided
we can find the required values of $\alpha$ and $\beta$ from data (i.e.
simulation). In the same figure we show fragmentation functions for lower $M$,
which demonstrates that the produced glueball spectrum becomes dominated by
finite-mass effects as $M\to 10m_{0}$, since small $x$-values become
kinematically forbidden and a significant fraction of events only produce two
glueballs with $x=1$. This again demonstrates the limitations of using
analytical approximations to estimate glueball distributions, and underlines
the importance of using a self-consistent Monte Carlo simulation like
GlueShower.

We also show the exclusive fragmentation functions of the lightest 10 glueball
stats for $M/m_{0}=100$ in Fig. 5 (right). For $T_{\text{had}}=T_{c}$, the
final states predominately consist of $0^{++}$ and $2^{++}$ glueballs,
followed by the next three heaviest states in roughly equal proportion.
Heavier state production is suppressed by an order of magnitude.

###  V.2 Plasma-like Hadronization

We now discuss some of the same GlueShower outputs for a basic plasma-like
shower case. The orange curve in Fig. 4 shows how multiplicity scales with
initial center-of-mass energy for $\Lambda_{\text{had}}/(2m_{0})=c=4$. Overall
we find that the plasma-like case leads to higher multiplicity events compared
to the jet-like case, but also that the inclusive multiplicity scales
differently in various energy regimes:

  * •

In the large energy limit $M\gg\Lambda_{\text{had}}$, we recover the same
multiplicity scaling as the jet-like case, but larger by a constant since each
produced plasma-ball evaporates into at least two glueballs. This does not
mean, however, that the plasma-like case converges to the jet-like case in the
high-energy limit: while the overall multiplicity scales similarly with
energy, the absolute multiplicity is higher in the plasma-like case, resulting
in lower characteristic energies for the final-state glueballs.131313It is
helpful to consider a SM QCD analogy: one could imagine a hypothetical
hadronization model that dominantly produces heavy $B$-mesons (analogous to
the plasma-balls), which then decay to lighter hadrons. Certainly, various
multiplicity and momentum scalings of this $B$-meson-shower would be similar
to the scalings of pion multiplicity and momenta in realistic SM
hadronization, but the different hadronization assumption would also introduce
fudamental differences in the final states no matter what the initial energy.

  * •

For $M\to\Lambda_{\text{had}}$ the multiplicity asymptotes to a value larger
than 2, in this case roughly 5. This is the regime in which the initial state
only splits into two singlet plasma balls, which then each evaporate to at
least two glueballs each.

  * •

At $M=\Lambda_{\text{had}}$ there is a discontinuity as the shower enters a
new regime in which there is insufficient center-of-mass energy to form two
separate plasma balls. For $M\leq\Lambda_{\text{had}}$, the initial di-gluon
production event is therefore taken to form just a single colour-singlet
plasma ball of mass $M$, which evaporates by glueball emission. Since this
leads to a larger $m_{\text{plasma}}$ just below $\Lambda_{\text{had}}$ than
just above, there is a small spike in produced glueball multiplicity, but this
is a reasonable threshold effect. In this regime, glueball production is
entirely SUEP-like, leading to a very different multiplicity scaling until the
absolute minimum of $N=2$ is reached for $M<3m_{0}$.

Inclusive and exclusive glueball fragmentation functions in the plasma-like
case are shown in Fig. 6. As demonstrated in the left plot, the analytical
functional form of the inclusive fragmentation function is still a fair
approximation in the high-energy large-$x$ regime, though significantly worse
than for the jet-like case, and for modest or low energies, the differences
are even more pronounced. On the right we show exclusive energy spectra for
the lightest 10 glueball species, which are produced with very similar
relative multiplicities as in the jet-like case. This is what we would expect,
as relative glueball multiplicities are determined by the hadronization
temperature which is kept constant between the cases considered in both Fig. 5
and Fig. 6. Additionally this shows that threshold effects in the jet-like
case that favour the production of light glueballs have minimal effect.

###  V.3 Defining Hadronization Benchmarks

Figure 7: Effect of $c=\Lambda_{\text{had}}/(2m_{0})$ on average glueball
multiplicity for exotic Higgs decays ($M=125\;\mathrm{GeV}$) into glueballs
with $m_{0}=10\;\mathrm{GeV}$. In the jet-like case, $\Lambda_{\text{had}}$ is
the hadronization scale at which the shower is terminated and gluons turned
into glueballs. In the plasma-like case, it is the twice the mass of colour-
singlet plasma balls produced in the shower, which then evaporate into
glueballs. In both cases we take $T_{\text{had}}/T_{c}$=1. The red dashed line
indicates $\Lambda_{\text{had}}=M$, resulting in the production of just two
glueballs in the jet-like case and a single plasma ball of mass $M$ in the
plasma-like case.

Having discussed how some simple observables behave at different energies in
both the jet-like and plasma-like cases, we now systematically examine their
dependence on two nuisance parameters of our simulation: the hadronization
scale $\Lambda_{\text{had}}=c\cdot 2m_{0}$ and the hadronization temperature
$T_{\text{had}}=d\cdot T_{c}$. This allows us to argue for a small set of
benchmark hadronization parameters that should span the range of physically
reasonable possible outcomes for glueball production, and which hence define
theory uncertainties for physical predictions.

Figure 8: Dependence of average glueball multiplicity and energy on $M/m_{0}$
for the 8 hadronization benchmarks with $N_{c}$ = 3. The range of predictions
can be interpreted as our current theoretical uncertainty on glueball
production: a factor of a few on average glueball energy and inclusive
multiplicity, and a factor of 10 on exclusive $0^{++}$ multiplicity.  Figure
9: $0^{++}$ glueball energy spectra for $M/m_{0}=10$ (left) and 100 (right)
and $N_{c}=3$, for our 8 hadronization benchmarks.

Figure 7 shows how the inclusive glueball multiplicity in both the jet- and
plasma-like cases depends on $\Lambda_{\text{had}}$ for
$T_{\text{had}}=T_{c}$. In the jet-like case, shown in blue, we see that as
$\Lambda_{\text{had}}$ increases the multiplicity decreases and asymptotes to
2. The plasma-like case (shown in orange) requires a minimum value of $c=2$ to
allow each plasma ball to emit at least two glueballs. At this limit, only the
lightest glueball can be produced. Naively we would expect multiplicity to
increase with $c$, because even as increasing $c$ logarithmically suppresses
the number of splittings in the shower, it linearly increases the number of
glueballs emitted by the larger plasma states. Somewhat surprisingly, from
$c=$ 2 to $\sim$ 3 the multiplicity decreases. This is due to the fact the
plasma still only decays to two glueballs most of the time across this range,
but now with greater access to the higher mass glueball states. Not until
$c\sim 4$ is the entire glueball spectrum sampled and produced in plasma ball
decays, and further increasing $c$ now does lead to an increase in glueballs
produced by the plasma. Finally, as discussed in the previous subsection, when
$\Lambda_{\text{had}}\to M$ only a single plasma ball is produced in the
event, resulting in fully SUEP-like behaviour.

While Fig. 7 depicts both hadronization interpretations across the full range
of $c$ values, each case is physically sensible in separate regimes. The jet-
like interpretation assumes the initial flux tube fragments directly into
final state glueballs, which is only sensible for values of $c$ close to 1.
Larger values force gluons with virtualities $\sqrt{t}\gg m_{0}\sim 6\Lambda$
into a single low-mass glueball, which seems very implausible given the a
priori reliable perturbative prediction for a much higher jet mass of that
branch of the shower. We therefore adopt $c=1,2$ as benchmark values for the
hadronization scale in the jet-like regime.

In contrast, the plasma-like case assumes fragmentation into high mass plasma
states, and small values $c\lesssim 4$ seem to imply that these plasma balls
have a large bias of evaporating only into the light glueball states. Sampling
the full range of plasma-like behaviour therefore motivates picking slightly
larger values of $c$. On the other hand, the plasma-like scenario already
defies the most reasonable physical expectation based on our understanding of
perturbative QCD and flux-tube dynamics, and a very large value of the plasma
ball mass pushes this scenario into even more implausible regimes. For
example, $c=6$ corresponds to ending the perturbative shower at $\sim
70\Lambda$, where $SU(N_{c})$ should be entirely perturbative. It should
therefore serve as a suitable ceiling for the possible range of plasma masses
produced under this exotic assumption for the non-perturbative behaviour of
$N_{f}=0$ QCD. We therefore adopt $c=4,6$ as benchmark values for the
hadronization scale in the plasma-like regime.

The hadronization temperature $T_{\text{had}}=d\cdot T_{c}$ is less
constrained by physical arguments and relatively unimportant for inclusive
observables across glueball species. However, $T_{\text{had}}$ dominantly
determines the exclusive predictions for each type of glueball. Given the
exponential dependence of both glueball momentum in the plasma-like case and
relative multiplicities in both jet- and plasma-like cases on
$T_{\text{had}}$, sampling $d=1,2$ should span a wide range of physically
plausible predictions.

In summary, careful phenomenological studies involving dark glueball
production should compute physical predictions for a variety of different
values of the nuisance parameters in GlueShower. For the jet-like case,
$c=1,2$ and $d=1,2$ should be simulated. For the plasma-like case, $c=4,6$ and
$d=1,2$. To be conservative, all 8 benchmark points should be used to define
the systematic error bar on predictions.

###  V.4 New predictions for Glueball Production

We now have in hand a physically motivated simulation of the glueball
production process, as well as 8 hadronization benchmarks which span the range
of physically reasonable possible outcomes. This allows us to make fully self-
consistent predictions for glueball production with accurate theoretical
uncertainties included for the first time.

Naively, we would expect collider signals of dark glueball production to be
most sensitive to the multiplicity and decay mode of the shortest-lived or
most visibly-decaying glueball state, for example the $0^{++}$ if decay
proceeds via the Higgs portal. Conversely, indirect detection of dark matter
annihilating into dark glueballs Curtin and Gemmell  will be affected by the
relative distributions and decays of all unstable glueball species. To give a
feeling for how each of these two types of studies might be affected by
theoretical uncertainties, we show in Fig. 8 how average multiplicity and
energy predictions for all glueballs inclusively and for the $0^{++}$
exclusively change across our range of possible hadronization benchmarks. We
take the variation across all benchmarks to indicate the theoretical
uncertainty for each observable.

We find that average glueball multiplicity, average glueball energy and
average $0^{++}$ energy have an uncertainty of about a factor of 3 across the
range of considered center-of-mass energies $M/m_{0}$. On the other hand, the
exclusive multiplicity of the $0^{++}$ state has a much larger spread of
possible predictions, roughly a factor of 10 across the hadronization
benchmarks. Slightly more can be said if one is willing to ascribe different
priors to the default jet-like versus the more exotic plasma-like
hadronization hypotheses. Jet-like showers produce lower multiplicity jets of
higher energy glueballs, while plasma-like showers produce fatter jets with
higher multiplicities of softer glueballs.

Figure 9 compares the $0^{++}$ energy spectra for different hadronization
benchmarks. At low center-of-mass energies, we find very significant
differences between the jet-like and plasma-like showers. While the plasma-
like case favours low energy glueballs, the jet-like case is roughly flat,
with dominant contributions by simple two-body glueball production. At high
center-of-mass energies, the spectral shapes become much more similar, with
low energy final states favoured, to varying degrees, across the benchmarks.

##  VI Conclusion

Hidden Valleys are an extremely popular hypothesis for BSM physics. They may
address fundamental mysteries like the hierarchy problem or the nature of dark
matter, and their signatures are targeted by many new LHC searches (e.g.
Alimena  _et al._ (2020); Aad  _et al._ (2013, 2020); Sirunyan  _et al._
(2019); Tumasyan  _et al._ (2021)) and proposed detectors Curtin  _et al._
(2019); Alpigiani  _et al._ (2020); Feng  _et al._ (2018); Gligorov  _et al._
(2018). However, the seemingly simple and minimal case of $N_{f}=0$ hidden QCD
has undergone very little quantitative study, despite its high theoretical
motivation within frameworks like Neutral Naturalness Chacko  _et al._
(2006a); Craig  _et al._ (2015); Burdman  _et al._ (2007); Barbieri  _et al._
(2005); Chacko  _et al._ (2006b); Cai  _et al._ (2009); Poland and Thaler
(2008); Cohen  _et al._ (2018); Cheng  _et al._ (2018). This can be traced
back to our lack of understanding of pure glue hadronization.

In this work, we show that significant progress can be made by combining a
perturbative pure glue parton shower with a self-consistent and physically
motivated parameterization of the unknown non-perturbative physics. This is in
large part enabled by the modest hierarchy between the glueball mass and the
confinement scale $m_{0}/\Lambda\sim 6$ in $SU(N_{c})$ theories. We make our
simulation code available as the public Python code GlueShower, the first
glueball generator for Hidden Valley theories.141414GlueShower download:
github.com/davidrcurtin/GlueShower This can be combined with other Monte Carlo
generators for dark gluon production and dark glueball decay Juknevich (2010);
Juknevich  _et al._ (2009) for a complete simulation of glueball production
within a given dark sector scenario.

Our approach incorporates what we consider to be an almost maximally wide
range of possible hadronization scenarios, to make sure no physically
reasonable possibility escapes our study. We argue that glueball production
proceeds largely similar to jet-like hadron production in SM QCD, but we also
consider the possibility that highly exotic non-perturbative physics of pure
Yang-Mills theories somehow favours the production of high-mass colour-singlet
gluon-plasma states, which evaporate via isotropic emission of thermal
glueballs in their restframe. In both jet-like and plasma-like modes, a range
of colour-singlet formation scales and hadronization temperatures can be
selected to generate a range of possible hadronization behaviours. We intend
GlueShower to be a starting point for $N_{f}=0$ dark sector studies, and point
out future improvements that could be implemented. A careful treatment of
colour confinement or a study of gluons in the rope-like picture Biro  _et
al._ (1984) could lead to a more sophisticated hadronization model, but the
current version represents a physically reasonable starting point for
phenomenological studies.

We define a set of 8 hadronization benchmark scenarios, 4 in each of the above
two modes, which we argue span the range of physically remotely reasonable
hadronization outcomes for glueball production. We then use GlueShower to make
some preliminary studies of glueball production observables, with theoretical
uncertainties obtained from the variation across the different hadronization
benchmarks. Uncertainties range from modest $\mathcal{O}(1)$ factors to orders
of magnitude, depending on the observable.

Accurately taking these uncertainties into account will be important for
designing and interpreting future searches and constraints. In particular, our
work could highlight which observables are more robust with respect to unknown
details of glueball hadronization, encouraging a variety of collider and
astrophysical probes to assemble a complete picture of the dark sector
fundamental parameters and nonperturbative dynamics. It is our hope that this
work enables many new studies and searches of dark sectors that were
previously intractable.

_Acknowledgements:_ We especially thank Matt Strassler for many insightful
conversations and discussions. We also thank Jared Barron and Benjamin Fuks
for helpful conversations. The research of DC and CG was supported in part by
a Discovery Grant from the Natural Sciences and Engineering Research Council
of Canada, the Canada Research Chair program, the Alfred P. Sloan Foundation,
and the Ontario Early Researcher Award. The work of CG was also supported by
the University of Toronto Connaught International Scholarship.

##  Appendix A Perturbative QCD shower

We provide a step-by-step description of how GlueShower generates perturbative
showers of gluons. This is a completely standard lowest-order perturbative
shower with angular ordering, implemented following the Pythia manual
Sjostrand  _et al._ (2006), but we include this summary here for completeness
and pedagogical purposes.

The two initial gluons are treated separately to ensure energy-momentum
conservation between their initial virtualities and energies. The shower is
then evolved iteratively for all the subsequent daughter dark gluons.

###  A.1 Initial Gluon Setup

  1. (1)

Find virtuality of gluon 1 ($t_{1}$) assuming gluon 2 is on-shell. This is
done by starting at the maximum allowed value, $(M-m_{\text{min}})^{2}$, and
evolving down using the Monte Carlo method described in Sec. III. There is
some ambiguity regarding what on-shell means when the final state glueball
species is still unspecified. For consistency across the shower, we define the
minimum mass as half the hadronization scale, $c\cdot m_{0}$.

  2. (2)

Once $t_{1}$ is found, if this allows the gluon to split, $z_{1}$ is
determined using the value of $t_{1}$.

  3. (3)

Steps 1 and 2 are repeated for gluon 2, to find $t_{2}$, and possibly $z_{2}$
if gluon 2 also splits.

  4. (4)

Using these values, a range of vetoes are checked before continuing with the
shower.

     * •

First check that $t_{1}+t_{2}<M$, if false then evolve the larger of the two
virtualities to a smaller value.

     * •

$z_{1,2}$ are found assuming that the other gluon was on-shell, but now they
have virtualities of $t_{1,2}$; check that $z_{1,2}$ remain within the updated
allowed range. If a gluon fails this check it virtuality is evolved to smaller
values.

     * •

Above vetoes are repeated until passed.

  5. (5)

We now have values for the virtuality and splitting energy of each initial
gluon, which are then used to evolve the daughters. Note that the 4-vectors of
the initial gluons are now entirely determined.

###  A.2 Evolving the Daughters

We now outline the steps applied for general gluon evolution at any point
along the shower, excluding the initial gluons. We label the produced daughter
gluons in this step 3 and 4.

  1. (1)

First, the coordinate system of this new splitting is established:

     * •

The $z$ axis is aligned with the direction of the parent gluon’s momentum.

     * •

A random angle is chosen in the $x-y$ plane for the perpendicular momentum of
the daughters.

  2. (2)

Find initial $t$ and $z$ guess for the daughters.

     * •

Initial energies are given by, $E_{3}=z_{1}*E_{1}$ and
$E_{4}=(1-z_{1})*E_{1}$.

     * •

The starting virtuality, $t_{\text{start},i}$ , for each daughter is given by
min($E_{i}^{2},m_{1}^{2})$.

     * •

Evolve the daughters down to find first guesses for $t_{3}$ and $t_{4}$. If
splitting is allowed use these values to find $z_{3}$ and $z_{4}$.

  3. (3)

Apply unconstrained/constrained evolution for each daughter. In unconstrained
evolution, the maximum possible kinematic $z$ range is given when assuming the
daughters have virtuality $m_{\text{min}}^{2}$. Thus, when the daughters are
evolved and their actual virtualities are found, $t_{3,4}$, the originally
assumed $z$ value may lie outside the newly determined range. The solution is
to redefine the $z$ value. This is done by taking the original four momenta of
the daughters and boosting them to their centre of mass frame, rescaling the
vectors to reflect their determined evolved virtuality, and then boosting them
back to the lab frame. This leads to a new $z$ value, matching the new
energies of the rescaled dark gluons, that lies within the required range and
is equivalent to using Eq. (10.11) given in Sjostrand  _et al._ (2006).

For constrained evolution, the final masses of the daughter gluons are
constrained by the requirement that $z$ remain in the originally determined
range. The default option enabled in the GlueShower code is unconstrained
evolution, as used in Pythia, and leads to a higher rate of splittings in the
shower.

  4. (4)

Impose angular ordering. This is the result of soft gluon coherence effects
which cause the opening angle between the daughter gluons to be smaller than
the opening angle of the parent gluon. This can be described as the supression
of wide angle gluon emission. Practically it is enforced by determining the
absolute maximum opening angle of the daughters, by assuming their daughters
have mass $m_{\text{min}}$, and if this angle is bigger than the parent’s
opening angle, the daughter is evolved to lower virtuality. Note that if a
gluon is evolved to satisfy angular ordering, step 4 is repeated to ensure
that the unconstrained evolution condition still holds.

  5. (5)

With the $t,z$ values of the daughter gluons determined, the steps of this
section are repeated iteratively until all dark gluons in the shower have
reached the hadronization scale.

## References

  * Strassler and Zurek (2007) M. J. Strassler and K. M. Zurek, Phys. Lett. B 651, 374 (2007), arXiv:hep-ph/0604261 .
  * Kang and Luty (2009) J. Kang and M. A. Luty, JHEP 11, 065 (2009), arXiv:0805.4642 [hep-ph] .
  * Bai and Schwaller (2014) Y. Bai and P. Schwaller, Phys. Rev. D 89, 063522 (2014), arXiv:1306.4676 [hep-ph] .
  * Renner and Schwaller (2018) S. Renner and P. Schwaller, JHEP 08, 052 (2018), arXiv:1803.08080 [hep-ph] .
  * Mies  _et al._ (2021) H. Mies, C. Scherb, and P. Schwaller, JHEP 04, 049 (2021), arXiv:2011.13990 [hep-ph] .
  * Holdom (1986) B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Patt and Wilczek (2006) B. Patt and F. Wilczek, (2006), arXiv:hep-ph/0605188 .
  * Falkowski  _et al._ (2009) A. Falkowski, J. Juknevich, and J. Shelton, (2009), arXiv:0908.1790 [hep-ph] .
  * Chacko  _et al._ (2006a) Z. Chacko, H.-S. Goh, and R. Harnik, Phys. Rev. Lett. 96, 231802 (2006a), arXiv:hep-ph/0506256 .
  * Craig  _et al._ (2015) N. Craig, A. Katz, M. Strassler, and R. Sundrum, JHEP 07, 105 (2015), arXiv:1501.05310 [hep-ph] .
  * Burdman  _et al._ (2007) G. Burdman, Z. Chacko, H.-S. Goh, and R. Harnik, JHEP 02, 009 (2007), arXiv:hep-ph/0609152 .
  * Barbieri  _et al._ (2005) R. Barbieri, T. Gregoire, and L. J. Hall, (2005), arXiv:hep-ph/0509242 .
  * Chacko  _et al._ (2006b) Z. Chacko, Y. Nomura, M. Papucci, and G. Perez, JHEP 01, 126 (2006b), arXiv:hep-ph/0510273 .
  * Cai  _et al._ (2009) H. Cai, H.-C. Cheng, and J. Terning, JHEP 05, 045 (2009), arXiv:0812.0843 [hep-ph] .
  * Poland and Thaler (2008) D. Poland and J. Thaler, JHEP 11, 083 (2008), arXiv:0808.1290 [hep-ph] .
  * Cohen  _et al._ (2018) T. Cohen, N. Craig, G. F. Giudice, and M. Mccullough, JHEP 05, 091 (2018), arXiv:1803.03647 [hep-ph] .
  * Cheng  _et al._ (2018) H.-C. Cheng, L. Li, E. Salvioni, and C. B. Verhaaren, JHEP 05, 057 (2018), arXiv:1803.03651 [hep-ph] .
  * Alimena  _et al._ (2020) J. Alimena _et al._ , J. Phys. G 47, 090501 (2020), arXiv:1903.04497 [hep-ex] .
  * Aad  _et al._ (2013) G. Aad _et al._ (ATLAS), JINST 8, P07015 (2013), arXiv:1305.2284 [hep-ex] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 80, 450 (2020), arXiv:1909.01246 [hep-ex] .
  * Strassler (2008) M. J. Strassler, (2008), arXiv:0801.0629 [hep-ph] .
  * Knapen  _et al._ (2017) S. Knapen, S. Pagan Griso, M. Papucci, and D. J. Robinson, JHEP 08, 076 (2017), arXiv:1612.00850 [hep-ph] .
  * Barron  _et al._ (2021) J. Barron, D. Curtin, G. Kasieczka, T. Plehn, and A. Spourdalakis, (2021), arXiv:2107.12379 [hep-ph] .
  * Knapen  _et al._ (2021) S. Knapen, J. Shelton, and D. Xu, Phys. Rev. D 103, 115013 (2021), arXiv:2103.01238 [hep-ph] .
  * Cohen  _et al._ (2015) T. Cohen, M. Lisanti, and H. K. Lou, Phys. Rev. Lett. 115, 171804 (2015), arXiv:1503.00009 [hep-ph] .
  * Cohen  _et al._ (2017) T. Cohen, M. Lisanti, H. K. Lou, and S. Mishra-Sharma, JHEP 11, 196 (2017), arXiv:1707.05326 [hep-ph] .
  * Cohen  _et al._ (2020) T. Cohen, J. Doss, and M. Freytsis, JHEP 09, 118 (2020), arXiv:2004.00631 [hep-ph] .
  * Schwaller  _et al._ (2015) P. Schwaller, D. Stolarski, and A. Weiler, JHEP 05, 059 (2015), arXiv:1502.05409 [hep-ph] .
  * Linthorne and Stolarski (2021) D. Linthorne and D. Stolarski, Phys. Rev. D 104, 035019 (2021), arXiv:2103.08620 [hep-ph] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), JHEP 02, 179 (2019), arXiv:1810.10069 [hep-ex] .
  * Tumasyan  _et al._ (2021) A. Tumasyan _et al._ (CMS), (2021), arXiv:2112.11125 [hep-ex] .
  * Morningstar and Peardon (1999) C. J. Morningstar and M. J. Peardon, Phys. Rev. D 60, 034509 (1999), arXiv:hep-lat/9901004 .
  * Lucini and Moraitis (2008) B. Lucini and G. Moraitis, Phys. Lett. B 668, 226 (2008), arXiv:0805.2913 [hep-lat] .
  * Teper (1998) M. J. Teper, journal (1998), arXiv:hep-th/9812187 .
  * Lucini  _et al._ (2010) B. Lucini, A. Rago, and E. Rinaldi, JHEP 08, 119 (2010), arXiv:1007.3879 [hep-lat] .
  * Athenodorou and Teper (2021) A. Athenodorou and M. Teper, journal (2021), arXiv:2106.00364 [hep-lat] .
  * Yamanaka  _et al._ (2021a) N. Yamanaka, A. Nakamura, and M. Wakayama, in  _38th International Symposium on Lattice Field Theory_ (2021) arXiv:2110.04521 [hep-lat] .
  * Juknevich (2010) J. E. Juknevich, JHEP 08, 121 (2010), arXiv:0911.5616 [hep-ph] .
  * Juknevich  _et al._ (2009) J. E. Juknevich, D. Melnikov, and M. J. Strassler, JHEP 07, 055 (2009), arXiv:0903.0883 [hep-ph] .
  * Faraggi and Pospelov (2002) A. E. Faraggi and M. Pospelov, Astropart. Phys. 16, 451 (2002), arXiv:hep-ph/0008223 .
  * Boddy  _et al._ (2014a) K. K. Boddy, J. L. Feng, M. Kaplinghat, and T. M. P. Tait, Phys. Rev. D 89, 115017 (2014a), arXiv:1402.3629 [hep-ph] .
  * Boddy  _et al._ (2014b) K. K. Boddy, J. L. Feng, M. Kaplinghat, Y. Shadmi, and T. M. P. Tait, Phys. Rev. D 90, 095016 (2014b), arXiv:1408.6532 [hep-ph] .
  * Garcia Garcia  _et al._ (2015) I. Garcia Garcia, R. Lasenby, and J. March-Russell, Phys. Rev. D 92, 055034 (2015), arXiv:1505.07109 [hep-ph] .
  * Soni and Zhang (2016) A. Soni and Y. Zhang, Phys. Rev. D 93, 115025 (2016), arXiv:1602.00714 [hep-ph] .
  * Soni  _et al._ (2017) A. Soni, H. Xiao, and Y. Zhang, Phys. Rev. D 96, 083514 (2017), arXiv:1704.02347 [hep-ph] .
  * Forestell  _et al._ (2017) L. Forestell, D. E. Morrissey, and K. Sigurdson, Phys. Rev. D 95, 015032 (2017), arXiv:1605.08048 [hep-ph] .
  * Forestell  _et al._ (2018) L. Forestell, D. E. Morrissey, and K. Sigurdson, Phys. Rev. D 97, 075029 (2018), arXiv:1710.06447 [hep-ph] .
  * Yamanaka  _et al._ (2021b) N. Yamanaka, H. Iida, A. Nakamura, and M. Wakayama, Phys. Lett. B 813, 136056 (2021b), arXiv:1910.01440 [hep-ph] .
  * Yamanaka  _et al._ (2020) N. Yamanaka, H. Iida, A. Nakamura, and M. Wakayama, Phys. Rev. D 102, 054507 (2020), arXiv:1910.07756 [hep-lat] .
  * Jo  _et al._ (2021) B. Jo, H. Kim, H. D. Kim, and C. S. Shin, Phys. Rev. D 103, 083528 (2021), arXiv:2010.10880 [hep-ph] .
  * Spergel and Steinhardt (2000) D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84, 3760 (2000), arXiv:astro-ph/9909386 .
  * Weinberg  _et al._ (2015) D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray, and A. H. G. Peter, Proc. Nat. Acad. Sci. 112, 12249 (2015), arXiv:1306.0913 [astro-ph.CO] .
  * Curtin and Verhaaren (2015) D. Curtin and C. B. Verhaaren, JHEP 12, 072 (2015), arXiv:1506.06141 [hep-ph] .
  * Chacko  _et al._ (2016) Z. Chacko, D. Curtin, and C. B. Verhaaren, Phys. Rev. D 94, 011504 (2016), arXiv:1512.05782 [hep-ph] .
  * Burdman and Lichtenstein (2018) G. Burdman and G. Lichtenstein, JHEP 08, 146 (2018), arXiv:1807.03801 [hep-ph] .
  * Andersson  _et al._ (1983) B. Andersson, G. Gustafson, G. Ingelman, and T. Sjöstrand, Physics Reports 97, 31 (1983).
  * Webber (1984) B. Webber, Nuclear Physics B 238, 492 (1984).
  * Bahr  _et al._ (2008) M. Bahr _et al._ , Eur. Phys. J. C 58, 639 (2008), arXiv:0803.0883 [hep-ph] .
  * Amati and Veneziano (1979) D. Amati and G. Veneziano, Physics Letters B 83, 87 (1979).
  * Carloni and Sjostrand (2010) L. Carloni and T. Sjostrand, JHEP 09, 105 (2010), arXiv:1006.2911 [hep-ph] .
  * Carloni  _et al._ (2011) L. Carloni, J. Rathsman, and T. Sjostrand, JHEP 04, 091 (2011), arXiv:1102.3795 [hep-ph] .
  * Bierlich  _et al._ (2022) C. Bierlich _et al._ , (2022), arXiv:2203.11601 [hep-ph] .
  * (63) D. Curtin and C. Gemmell, work in progress .
  * Prosperi  _et al._ (2007) G. M. Prosperi, M. Raciti, and C. Simolo, Prog. Part. Nucl. Phys. 58, 387 (2007), arXiv:hep-ph/0607209 .
  * Batell  _et al._ (2020) B. Batell, W. Hu, and C. B. Verhaaren, JHEP 08, 009 (2020), arXiv:2004.10761 [hep-ph] .
  * Kilic  _et al._ (2021) C. Kilic, C. B. Verhaaren, and T. Youn, Phys. Rev. D 104, 116018 (2021), arXiv:2109.03248 [hep-ph] .
  * Hou  _et al._ (2001) W.-S. Hou, C.-S. Luo, and G.-G. Wong, Phys. Rev. D 64, 014028 (2001), arXiv:hep-ph/0101146 .
  * Ishii  _et al._ (2001) N. Ishii, H. Suganuma, and H. Matsufuru, in  _Workshop on Lepton Scattering, Hadrons and QCD_ (2001) pp. 252–259, arXiv:hep-lat/0106004 .
  * Loan and Ying (2006) M. Loan and Y. Ying, Prog. Theor. Phys. 116, 169 (2006), arXiv:hep-lat/0603030 .
  * Ellis  _et al._ (2011) R. K. Ellis, W. J. Stirling, and B. R. Webber,  _QCD and collider physics_ , Vol. 8 (Cambridge University Press, 2011).
  * Sjostrand  _et al._ (2006) T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 .
  * Marchesini and Webber (1984) G. Marchesini and B. R. Webber, Nucl. Phys. B 238, 1 (1984).
  * ’t Hooft (1974) G. ’t Hooft, Nucl. Phys. B 72, 461 (1974).
  * Bulava  _et al._ (2019) J. Bulava, B. Hörz, F. Knechtli, V. Koch, G. Moir, C. Morningstar, and M. Peardon, Phys. Lett. B 793, 493 (2019), arXiv:1902.04006 [hep-lat] .
  * Minkowski and Ochs (2000) P. Minkowski and W. Ochs, Phys. Lett. B 485, 139 (2000), arXiv:hep-ph/0003125 .
  * Field and Feynman (1978) R. Field and R. Feynman, Nuclear Physics B 136, 1 (1978).
  * Schwinger (1951) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Casher  _et al._ (1979) A. Casher, H. Neuberger, and S. Nussinov, Phys. Rev. D 20, 179 (1979).
  * Manes (2002) J. L. Manes, Nucl. Phys. B 621, 37 (2002), arXiv:hep-th/0109196 .
  * Blanchard  _et al._ (2004) P. Blanchard, S. Fortunato, and H. Satz, Eur. Phys. J. C 34, 361 (2004), arXiv:hep-ph/0401103 .
  * Noronha-Hostler  _et al._ (2010) J. Noronha-Hostler, J. Noronha, and C. Greiner, J. Phys. G 37, 094062 (2010), arXiv:1001.2610 [nucl-th] .
  * Petreczky (2012) P. Petreczky, J. Phys. G 39, 093002 (2012), arXiv:1203.5320 [hep-lat] .
  * Lucini  _et al._ (2004) B. Lucini, M. Teper, and U. Wenger, JHEP 01, 061 (2004), arXiv:hep-lat/0307017 .
  * Boyd  _et al._ (1996) G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier, and B. Petersson, Nucl. Phys. B 469, 419 (1996), arXiv:hep-lat/9602007 .
  * Lucini  _et al._ (2005) B. Lucini, M. Teper, and U. Wenger, JHEP 02, 033 (2005), arXiv:hep-lat/0502003 .
  * Lucini  _et al._ (2012) B. Lucini, A. Rago, and E. Rinaldi, Phys. Lett. B 712, 279 (2012), arXiv:1202.6684 [hep-lat] .
  * Biro  _et al._ (1984) T. S. Biro, H. B. Nielsen, and J. Knoll, Nucl. Phys. B 245, 449 (1984).
  * Bierlich (2018) C. Bierlich, EPJ Web Conf. 171, 14003 (2018), arXiv:1710.04464 [nucl-th] .
  * Swisdak (2013) M. Swisdak, Phys. Plasmas 20, 062110 (2013), arXiv:1305.5226 [physics.plasm-ph] .
  * Colangelo and Nason (1992) G. Colangelo and P. Nason, Phys. Lett. B 285, 167 (1992).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * de Florian  _et al._ (2017) D. de Florian, M. Epele, R. J. Hernandez-Pinto, R. Sassot, and M. Stratmann, Phys. Rev. D 95, 094019 (2017), arXiv:1702.06353 [hep-ph] .
  * Salas-Chavira  _et al._ (2021) A. Salas-Chavira, Z. Fan, and H.-W. Lin, (2021), arXiv:2112.03124 [hep-lat] .
  * Curtin  _et al._ (2019) D. Curtin _et al._ , Rept. Prog. Phys. 82, 116201 (2019), arXiv:1806.07396 [hep-ph] .
  * Alpigiani  _et al._ (2020) C. Alpigiani _et al._ (MATHUSLA), (2020), arXiv:2009.01693 [physics.ins-det] .
  * Feng  _et al._ (2018) J. L. Feng, I. Galon, F. Kling, and S. Trojanowski, Phys. Rev. D 97, 035001 (2018), arXiv:1708.09389 [hep-ph] .
  * Gligorov  _et al._ (2018) V. V. Gligorov, S. Knapen, M. Papucci, and D. J. Robinson, Phys. Rev. D 97, 015023 (2018), arXiv:1708.09395 [hep-ph] .