Dark Sector Glueballs at the LHC

We study confining dark sectors where the lightest hadrons are glueballs. Such models can provide viable dark matter candidates and appear in some neutral naturalness scenarios. In this work, we introduce a new phenomenological model of dark glueball hadronization inspired by the Lund string model. This enables us to make realistic predictions for dark glueball phenomenology at the LHC for the first time. Our model reproduces the expected thermal distribution of hadron species as an emergent consequence of hadronization dynamics. The ability to predict the production of glueball states heavier than the lightest species significantly expands the reach of long-lived glueball searches in MATHUSLA compared to previous simplified estimates. We also characterize regions of parameter space where emerging and/or semivisible jets could arise from pure-glue dark sectors, thereby providing new benchmark models that motivate searches for these signatures.

Previous collider studies have been limited by the absence of dedicated simulation tools for dark glueball production.Pythia's Hidden Valley module [79,80] is the current state-of-the-art in simulating strongly-coupled dark sector hadronization, but it does not accommodate the qualitatively different pure-glue (N f = 0) case.Recently, the Python module GlueShower [81] was created to address this gap.This package implemented a perturbative gluon shower and exploited only energy conservation and the large pure-glue mass-gap to parameterize the unknown details of glueball hadronization.While this enabled the first quantitative studies including the effects of the shower and multiple glueball species [71], the lack of any realistic hadronization dynamics resulted in very large uncertainties for exclusive quantities, like the production of specific glueball species, that can strongly influence collider signals.
In this work, we present a more sophisticated glueball hadronization model, based on applying Lund string dynamics [82] to the pure-glue regime, and implement it in a customized version of Pythia 8 [83].This enables more theoretically robust collider studies with improved uncertainties compared to the earlier GlueShower approach. 1e find that our hadronization algorithm dynamically realizes certain theoretically expected features, such as the thermal distribution of produced glueball species.These results suggest that the most important pure-glue hadronization dynamics may be captured by our approach.
We apply this implementation of dark glueball hadronization to classify the phenomenology across the parameter space of two specific models, determining which regions could possibly be probed by future emerging or semivisible jet searches.We consider applications to long-lived particle (LLP) searches and semivisible/emerging jet searches.In particular, we update the predicted signal at the proposed MATH-USLA LLP detector [84][85][86] of dark glueballs produced in exotic Higgs decays in Sec.4.1.Compared to earlier estimates based on two-body Higgs decays [44,87], our dark shower and hadronization dynamics lead to a dramatically expanded reach estimate, as the production of various glueball species with different lifetimes generates signals in different parts of parameter space.
Semivisible jets and emerging jets have become targets for LHC searches, see [19,20] and [22,23] respectively.This motivates developing consistent benchmark models that can serve as reference points for designing experimental analyses.In Sec.4.2, we show how glueball production can also realize both signatures and elucidate some properties of the resulting signals.In particular, the parameter r inv (which characterizes the collider-stable component of semivisible jets) can be predicted as a distribution, and we provide a number of examples in the results presented below.
The rest of this paper is structured as follows.We introduce the hadronization model in Sec. 2, first reviewing the Lund string model in Sec.2.1 and then showing our modifications of this approach for pure-glue hadronization in Sec.2.2 (further discussion is provided in App.A).We demonstrate that the expected thermal distribution of glueball species dynamically emerges in Sec.2.3.We provide benchmark values of the hadronization model parameters in Sec.2.4, which span a range of outputs representative of theoretical uncertainties from unknown hadronization details.In Sec. 3, we introduce glueball production and decay mechanisms through the Higgs portal that are relevant for collider signals.We then apply our glueball production and decay simulations to two phenomenology studies: glueball LLP decays in the proposed MATHUSLA experiment in Sec.4.1 and semivisible/emerging jet signatures in Sec.4.2.We conclude in Sec. 5.In App.B, we elaborate on the glueball species and momentum variations from our suggested benchmark parameters.A few additional distributions of interest for the semivisible jet scenario are provided in App. C.

Glueball Hadronization
The non-perturbative nature of hadronization necessitates introducing phenomenological models to make predictions for the collider signatures of confining sectors.This approach has a long history going back to Field and Feynman's "independent fragmentation" [88], where the model is simply that individual quarks fragment into the mesons making up a jet.The modern Monte Carlo event generator Pythia uses the "Lund string model" [82] (described below), while its contemporaries Herwig [89] and Sherpa [90] each use versions of "cluster hadronization" [91,92], where color-singlet clusters of partons that are close together in phase space decay into hadrons.Hadronization models introduce phenomenological nuisance parameters.
For SM QCD, one can perform elaborate tunes to data to constrain these parameters.Of course, we do not have the luxury of data to fit these parameters for a dark sector.Fortunately, as we show in Sec. 4, the observables of interest here are only moderately sensitive to the nuisance parameters.Developing theoretical predictions with error estimates and search strategies for dark sector jets that are less sensitive to nuisance parameters is an area of active interest, see e.g.[37,93].
The rest of this section is devoted to a description of string hadronization.First, we review the Lund string model used in Pythia 8.This provides the context for the description of our dark glueball hadronization model that follows.In all of our analyses, we consider the case where the number of colors N c is 3, but generalization to other values of N c is straightforward (given the glueball mass spectrum and its scaling with the confinement scale from lattice studies).

Lund String Model for Mesons and Baryons
We start by briefly summarizing the Lund string model of hadronization as implemented in Pythia 8.A more detailed explanation can be found in [83] with additional context in [82].Before hadronization, partons undergo a perturbative shower, iteratively splitting until the characteristic energy scale (transverse momentum p T relative to the parent parton in Pythia's implementation) reaches an IR cutoff p T min , which parameterizes the scale where the shower approaches strong coupling and must be matched onto the hadronization model.The shower is executed in the leading color (N c → ∞) approximation, such that each (anti-)quark has a unique (anti-)color label, and each gluon has unique color and anti-color labels.Therefore, for each color label, there is exactly one parton with the compensating anti-color label at each step of the shower.
Partons are grouped into color singlets by "Lund strings."These are simple representations of flux tubes dictating the flow of color charge.Quarks and antiquarks live at the ends of strings, while gluons are represented as kinks in the strings.Strings are therefore comprised of "string pieces," the segments of the string that connect individual quarks and gluons.Each string piece has momentum where p µ 1 and p µ 2 are the momenta of the partons connected by the string piece.This momentum defines a string piece mass m piece via p 2 piece = m 2 piece .Each (anti-)quark effectively donates all its momentum to the string piece that terminates on it, and each gluon donates half its momentum to each of the two string pieces that are connected to it.
In order to account for the fact that QCD has a finite number of colors, Pythia implements a procedure called "color reconnection" between the end of the perturbative shower and the formation of hadrons.Out of the few different options Pythia 8 offers for color reconnection, the so-called "QCD-based" version is the best-motivated   Color and anti-color labels are represented by filled and empty (semi-)circles, respectively.Quarks q and antiquarks q form the ends of strings, while gluons g are kinks in the string.String pieces span between the color/anti-color pairs of individual partons and are labeled with displayed colors (colors have numerical labels so that gray-scale versions contain the same information).Each parton is displayed with a fixed location in an abstract color-connection space where the displayed lengths of the string pieces correspond to the string-length λ. Between 1b and 1c, the pairs of blue (1) and green (3) connections are swapped.
for our regime without beam remnants or light quarks. 2Here, color/anti-color pairs are randomly reassigned a new label, which permits more than one possible grouping of partons into color singlets. 3Color reconnection seeks to minimize a Lorentzinvariant effective free energy λ called the string-length where m ref is the mass of some reference hadron.The minimization of λ is performed by swapping which color end of a color/anti-color pair is connected to which anticolor end.We provide a sketch of this procedure in Fig. 1 for a simple system with a few partons.The approach taken in the Lund string hadronization model is to consider these 2 Part of the purpose of color reconnection is to treat states at the end of showers and beam remnants consistently, so Pythia's default implementation was formulated with beam remnants specifically in mind.
3 More precisely, there are nine possible reassignments (with the restriction that gluons cannot be reassigned to have the same color and anti-color) to reflect the 1/9 probability of an SU (3) fundamental being able to form a singlet with an SU (3) anti-fundamental.More detail can be found in [83].
color connections between partons as oscillating classical strings that break up into hadrons.A string with two ends connecting a q q pair is called a "yo-yo mode" and is identified with a meson. 4When a hadron fragments off of a string, it takes away some random fraction z of the string's light-cone momentum p ± = E ± p z (with the z-axis being a preferred direction which Pythia takes to be the string axis), which has convenient Lorentz transformation properties.This z is sampled from the Lund Symmetric Fragmentation Function (LSFF) where a and b are phenomenological parameters, and the hadron's transverse mass m ⊥ appears due to quark tunneling effects explained in [83].Given that Lund fragmentation assumes strings break into pieces ending on quarks, this approach must be modified to accommodate a pure-glue sector.

String Model for Glueball Hadronization
In this section, we provide a qualitative discussion of our glueball hadronization algorithm for N c = 3, which can be easily generalized to other numbers of colors.Our focus is on explaining the role of the adjustable nuisance parameters that parameterize incalculable effects.A more detailed description with further justifications of the choices made here are is given in App. A. We modified Pythia 8 to simulate the branching of a Lorentz-singlet dark gluon pair (produced e.g. in the decay of a heavy scalar) via a leading color, p Tordered, pure-dark-glue parton shower.The shower cutoff scale is parameterized as p T min = c Λ D , where c is an O(1) nuisance parameter, and Λ D is the dark sector confinement scale. 6This cutoff scale parametrizes the onset of strongly-coupled dynamics, where the dark sector coupling α D becomes non-perturbatively large.In practice, the dimensionful scale that determines all of the scales in the pure-glue theory is the lightest glueball mass m 0 , as further described below.
Lattice studies have provided us with many glueball properties in pure SU (3), and in some cases for other values of N c [95][96][97].There are twelve species that are stable in the absence of external couplings, each with their own set of J P C quantum numbers [98].The lightest state is the 0 ++ with mass m 0 .Each heavy state's mass is a multiple of m 0 between ∼ 1.4 and ∼ 2.8, which we take from [96].The masses and spins of this glueball spectrum provide the inputs we will need to select the species of glueballs that are emitted during fragmentation, which will be described below.Lattice results also allow us to specify the boundary condition in the dark strong coupling's renormalization group evolution given a choice of m 0 , since m 0 = 6.28ΛD in pure SU (3) [97].In this way, we derive Λ D from the physical scale m 0 .We set the default value of c so that α D evaluated at the default shower cutoff scale is 1.
Non-default values of c change p T min without affecting the running of α D .Following the termination of the perturbative shower, we implement a version of QCD-based color reconnection.As visualized in Fig. 2, the only color-singlet Lund string topology in the pure-glue model is a closed ring.As in the SM, we swap color connections to minimize the string-length in Eq. (2.2), using the lightest glueball mass m 0 as m ref . 7This leaves us with color-singlet rings of string pieces that will fragment into glueballs.Both color reconnection and our hadronization algorithm are phenomenological models for flux tube rings twisting until they cross themselves and pinch into smaller rings.Without quarks, the strings are unable to break, so this pinching action is the only way the rings can divide into units with smaller invariant mass.
The intuitive picture of glueball rings pinching off the color-singlet flux tube ring in order to most rapidly decrease the total string-length inspires the glueball hadronization algorithm depicted in Fig. 3. First, we select string pieces with sufficient total invariant mass to be converted (or "fragmented") into a glueball.We then determine the species of the glueball by randomly selecting from among the species with mass less than that of the fragmenting string pieces, with probabilities weighted by the number of spin degrees of freedom.The emitted glueball's direction is along the total momentum of the fragmenting string pieces in the ring's rest frame, and its light-cone momentum is a random fraction z of the light-cone momentum of the fragmenting strings.We study the effect of sampling z from one of two fragmentation functions, the first being the LSFF in Eq. (2.3) with m ⊥ replaced by the glueball  mass m G , and the second being a beta distribution where α and k β are nuisance parameters.Having specified the glueball's species, direction, and light-cone momentum, its four momentum is fully determined by the on-shell condition.Whatever momentum is left over from the fragmenting string pieces is distributed equally among two new string pieces, and the resulting new ring can emit the next glueball.If emitting a glueball with the selected momentum would result in a new ring that is kinematically forbidden from further fragmentation, the ring instead fragments into two glueballs with the second glueball's species randomly sampled as previously discussed.In summary, our model's nuisance parameters are the shower cutoff factor c, the fragmentation function, and the chosen function's two shape parameters.

Emergence of Thermally Distributed Production Rates
By analogy to the hadron spectra produced in SM QCD fragmentation [99], a motivated expectation for the relative distribution of different glueball species is that they approximately follow a Maxwell-Boltzmann distribution [14]: where P J is the relative rate of producing the species with mass m J and spin J, and T had is some "hadronization temperature."In [14], which analyzed glueball production in a dark matter indirect detection context without separating out the perturbative shower, T had was taken to be related to the center-of-mass energy of the initial hard process.For high enough initial energy, this would result in a power law distribution favoring heavier hadron masses, a behavior quite different from what we are used to in SM QCD.Instead, we would expect that the hadronization temperature is fairly unrelated to the high scale at which the original gluons are produced (provided this high scale is sufficiently greater than any kinematic thresholds), but rather is set by an intrinsic feature of the confining theory.In particular, for N f = 0, N c = 3, the critical or Hagedorn temperature of the QCD phase transition [100][101][102] is T c ≃ 1.2 Λ D [97,103].This motivated the approach taken by GlueShower [81], which explicitly imposed this thermal distribution with T had = dT c for d ∼ O(1).Remarkably, we will show that our hadronization algorithm produces an approximately thermal multiplicity distribution of glueball species with very little dependence on the choice of fragmentation function or its parameters.Furthermore, the hadronization temperature is of the theoretically expected size, T had ≃ Λ D , increasing slowly with increasing shower cutoff scale for reasonable values of c ∼ O(1).
Our algorithm makes the minimal assumption that a local set of string pieces fragmenting into glueballs has no preference amongst the kinematically accessible glueball species, beyond the 2J + 1 spin multiplicity factor.The overall distribution of glueball species must therefore emerge from the kinematics of color-singlet rings.A color-singlet flux tube ring made of only soft string pieces that are close together in momentum space will predominantly produce only the lighter species, since each additional selected string piece will only modestly increase the total invariant mass available for fragmentation until the m 0 threshold is reached.On the other hand, fragmenting string pieces that are heavy compared to m 0 , as well as combinations of color-connected string pieces that are far-separated in momentum space, will produce heavy and light glueballs without preference.The combination of these effects results in a net suppression of the heavy species.This idea of "closeness" in momentum space (as determined by the invariant mass of sums of string piece momenta) elicits an intuitive geometric picture.If we imagine the color-singlet rings as polygons (as in Figs. 1 to 3) whose side lengths are determined by the string-length in Eq. (2.2), but whose angles randomly fluctuate, then combinations of string pieces that are "larger" in this sense have a greater propensity to cross each other and fragment off heavier glueballs.We would expect such a system to fragment in the order that most rapidly decreases its perimeter, which inspires our choice to begin fragmentation by selecting string pieces with the largest string-length. 8This intuitive picture may serve as a good analogy for the dynamics of closed flux tube rings.
This illustrates qualitatively how our hadronization algorithm provides a plausible model of glueball fragmentation, but there is no a priori reason to expect it to quantitatively produce an approximately thermal distribution.Furthermore, the above arguments suggest a significant dependence on the shower cutoff scale p T min = c Λ D , with higher values of c producing fewer gluon splittings and therefore fewer string pieces that each have higher mass, resulting in overproduction of heavy glueball states.Indeed, we observe a modest increase of the corresponding hadronization temperature with increasing c.
We now quantitatively demonstrate how this thermal species distribution emerges, and investigate the extent to which the glueball multiplicity distribution depends on the nuisance parameters of our hadronization model.Recall that these parameters are the shower cutoff scale set by c, the choice of fragmentation function between Eq. (2.3) or Eq.(2.4), and the two shape parameters for each function.Here, we focus on the fragmentation function and therefore set c to its default of 1.8.The discussion will not change for other O(1) values, and when defining hadronization benchmarks in the next section, we will include different choices for c.
The physical interpretation of the numerical values of the fragmentation function parameters is obscure, and there is no obvious correspondence between the parameters of the LSFF and that of the beta distribution.Therefore, we specify the mean µ z and standard deviation σ z of the probability distributions for the 0 ++ species, which fix values for the fragmentation function parameters and are easy to interpret.Another advantage of specifying the mean and standard deviation is that for each possible mean of a finitely-supported probability distribution, there is a maximum possible standard deviation.Thus, the space of all possible fragmentation function parameters is bounded when expressed this way.
In Fig. 4, we present the results of fitting the distributions of species to Eq. (2.5) for many points in the µ z -σ z plane for the fragmentation function of the 0 ++ species.We find the best fit T had by maximizing where y i is the total fraction of the i th species produced by the Monte Carlo, P J (x i ) is the thermal distribution in Eq. (2.5) evaluated at the mass x i of the i th species, and ⟨y⟩ is mean of the y i 's.We chose R 2 to quantify the quality of the fit, rather than χ 2 , because we are interested in the infinite Monte Carlo statistics limit.Using χ 2 would therefore artificially assign greater weight to the smallest species fractions.
The point with the best fit was in the plane of the beta distribution at µ z = 0.6, σ z = 0.008 with R 2 = 0.93 and T had = 0.99Λ D .We show results for shower center-of-  mass energy 125 GeV and glueball mass m 0 = 10 GeV, but the results are very similar for center of mass energy of 1 TeV.It is encouraging that the best-fit hadronization temperature lies close to the theoretical expectation, and that both T had and the high quality of fit depend only very little on the choice of fragmentation function, its parameters, or the center of mass energy.This suggests that our hadronization model may represent a good analogy for the true non-perturbative fragmentation dynamics of crossing color strings.
Fig. 5 shows a representative example of these production rates as a function of species, with the corresponding best fit using the default nuisance parameters described in Sec.2.4.Overall, the agreement between the thermal expectation and our model is fairly good.There is, however, a noticeable and consistent overproduction of very heavy glueball states compared to the Boltzmann expectation.
As an instructive comparison, Fig. 6 shows production rates of different SM hadron species produced by Pythia's default hadronization tune.The distribution organizes itself by the heaviest quark flavor in each species, and each of these groups Figure 5: Production rates P J of glueball species in ascending order of mass produced using the default parameters described in Sec.2.4 and the corresponding best fit to the thermal model in Eq. (2.5).The upper and lower plots have identical information, with the upper plot on a logarithmic scale and the lower plot on a linear scale.
appears to follow its own scaling relation as a function of mass.We fitted each group to its own Boltzmann distribution, generally finding good agreement between the fits and Monte Carlo, up to a few outliers.In fact, the heaviest SM hadrons appear to exhibit a mild enhancement compared to the thermal fit, which is consistent with the output of our glueball hadronization algorithm.These fits demonstrate that glueball species production from our algorithm should be approximately Boltzmann-like, but some deviations are to be expected, and the parameters should not necessarily be tuned solely to achieve an optimal Boltzmann fit.Rather, we are encouraged that our model's parameters only mildly impact the fit quality, and we set benchmarks in Sec.2.4 to capture the most extreme possible variations of the model's output.m (GeV)

Benchmark Parameters
In this section, we suggest three benchmarks for setting the shower cutoff and fragmentation function parameters for collider studies.Our analysis of thermal production rates in Sec.2.3 does not strongly favor a particular region in fragmentation function parameter space.Therefore, in addition to a well-motivated default choice, we define two bracketing variations that capture the different plausible outcomes of glueball hadronization.These are "soft" and "hard" scenarios, where the glueballs tend to have smaller and larger momentum in the dark shower rest frame, respectively.First, we investigate how varying the hadronization scale parameter c changes the glueball hardness.As the hadronization scale decreases, the gluons in the perturbative shower branch more often, and the Lund strings are softer, so the glueballs emitted by combining those strings are also softer.Thus, benchmarks for softer scenarios should correspond to smaller values of c.As for the fragmentation function parameters, it is simplest to interpret their effect in the µ z -σ z plane as in Sec.2.3, since parameters corresponding to a larger µ z tend to take more momentum from the fragmenting string pieces, resulting in harder glueballs.It is not obvious a priori which of the two fragmentation functions Eq. (2.3) and Eq.(2.4) would result in harder glueballs, but simulations show that the LSFF leads to harder kinematics by a small margin.See further discussion of glueball hardness in App.B.
With these qualitative effects in mind, we make some concrete suggestions for benchmarks.The default value of c is taken to be 1.8, since this sets our p T min at the scale where α D = 1.Pythia's default settings provide useful guidance in the variation of c because the SM α S evaluated at the default p T min for final state radiation is ≃ 1.6.The value of c for which our α D satisfies the same condition is 1.4, which we therefore adopt as our soft benchmark.We choose our hard benchmark c to be 2.1 so that the default is equidistant from the two variations.For the default fragmentation function, we choose the LSFF because this is what is used in the Lund string model.We choose the default function parameters to satisfy µ z = 0.5 and σ z = 0.3 since this is about as close to a uniform distribution as possible with this fragmentation function, and we do not have any reason to favor large or small z.For the soft and hard scenarios, we want to capture the extremes of the possible variations, so we choose µ z = 0.1 and 0.9, respectively, with σ z = 0.01.These numerical values, and their corresponding fragmentation function parameters, are summarized in Table 1.

Total Multiplicity Distributions
In addition to the relative production rates of different glueball species in Sec.2.3, we can compare distributions of the total number of glueballs per event N from our algorithm to analytical predictions for QCD with zero flavors.In the pure-glue limit with small hadron masses, the average number of hadrons is expected to scale with center of mass energy E CM as [104] ⟨N where C A is the quadratic Casimir factor for the adjoint representation, which is N c for SU (N c ). Fig. 7 shows a comparison of ⟨N ⟩ at various values of E CM between our algorithm with each set of benchmark parameters and the QCD prediction, with the normalization of Eq. (2.7) fixed by matching to the Monte Carlo at the largest E CM .The analytic prediction is slightly larger than our model at low energies, except where E CM is near 2m 0 and the prediction falls below the kinematic threshold of two hadrons.This is an expected consequence of the finite hadron masses, since for smaller E CM and parameters that tend to produce more glueballs, a larger portion of the energy is taken up by glueball masses.The reproduction of this standard result is a useful check on the validity of our algorithm.Notably, our algorithm tends to produce a greater number of glueballs than GlueShower, which generated ⟨N ⟩ ∼ 7 at E CM = 100 m 0 [81].Therefore, our more physically-motivated approach predicts greater discovery potential.

Dark Glueball Decay via Higgs Portal
We now have a concrete numerical method to simulate the production of dark glueballs.In order to connect with phenomenology, we need to specify the portal between the dark sector and the SM that will determine how the dark gluons are produced in a hard interaction and how the dark glueballs decay.We consider a pure-glue dark sector that couples to the SM via the Higgs portal, since this is the lowest-dimension portal that can connect the pure-glue sector to the SM.We compute the resulting glueball lifetimes and show them in the parameter space of two neutral naturalness scenarios.

Dark Glueball Decay Widths
In this section, we briefly summarize the pertinent results of [105] used in our estimation of dark sector glueball lifetimes and decay branching ratios.Strongly-coupled Glueball Mass (m 0 ) Higgs Portal dark sectors that include heavy fermions coupling to the Higgs give rise to the effective dimension-6 Higgs portal operator9 where H is the SM Higgs doublet, G µν (D) is the dark gluon field strength, α D is the dark sector strong coupling, M is the mass scale of the dark sector fermions, and y is an effective coupling that is determined by a model-dependent combination of the dark sector fermion Yukawa couplings with the Higgs (see [105] for explicit expressions).This operator can mediate both dark gluon production at the LHC and subsequent glueball decay to the SM.
The decay channels for each of the twelve glueballs are summarized in Table 2.The 0 ++ species decays into SM states ξ by mixing with the Higgs boson, 0 ++ → h * → ξξ, with the decay width where m h is the Higgs mass, Γ SM h→ξξ (m 0 ) is the decay width for a Higgs-like scalar of mass m 0 , which we calculate using HDECAY [106], and F 0 ++ is the non-perturbative decay constant with mass dimension 3. The heavier species (with the exceptions of the stable 0 −+ and 1 +− ) decay into lighter glueballs via emission of an off-shell Higgs, J → J ′ + h * (→ ξξ).The decay width for a glueball with spin J to a lighter glueball with spin J ′ and the SM is given by where g(x, y; z) = (1−x/z −y/z) 2 −4xy/z 2 , |M J,J ′ | is the averaged non-perturbative transition matrix element, and Γ (i) J,J ′ are dimensionless functions of the glueball masses that depend on the angular momentum transfer associated with each transition, which can be found in [105].The mass splitting of the glueball states can be a few GeV at small m 0 , where perturbative SM QCD breaks down, so we use Γ SM h→ξξ (m 12 ) values calculated from chiral perturbation theory [107] for this range.
The glueball decay widths depend on the theory parameters m 0 and M/y, as well as the non-perturbative decay constants and transition matrix elements.The matrix elements corresponding to the decays of the lightest glueballs have been calculated on the lattice [96], e.g.4π α D F 0 ++ = 2.3m 3 0 , which we use in this work.However, the matrix elements for the heavier states have not been computed.We use dimensional analysis to approximate the remaining heavy species' transition elements up to dimensionless prefactors, thus obtaining the correct scaling with m 0 .We set |M J,J ′ | = m 3 0 for our decay widths, vary each matrix element independently by a factor of two, and marginalize over the variation as part of our hadronization uncertainty.One could in principle also incorporate the dimension-8 operators listed in [105], which render the 0 −+ and 1 +− unstable.However, corrections due to these operators are suppressed by multiple orders of magnitude in the parameter space of interest, and these species only make up a few percent of produced glueballs, so we do not include these decays in our study.

Neutral Naturalness Models
The results of Sec.3.1 can be readily mapped onto parameters in neutral naturalness models because they generate the Higgs portal operator in Eq. (3.1) [44].For the Fraternal Twin Higgs (FTH) [45], where m tFTH is the twin top mass, and θ = tan −1 (m t /m tFTH ).For Folded Supersymmetry (FSUSY) [41], where m t is the SM top quark mass, m tFSUSY is the folded stop mass, and v is the SM Higgs vacuum expectation value.This mapping allows us to predict lifetimes and branching ratios for each glueball species given a choice of m 0 and the mass of the FSUSY or FTH top quark partner.Fig. 8 shows some representative examples.The plots show that 0 ++ state is the shortest lived state for each point in parameter space because it mixes directly with the Higgs.The next few heavy species that can decay via the dimension-6 operator (2 ++ , 2 −+ , 3 +− ) have much longer lifetimes due to only being able to radiate an off-shell Higgs.The remaining heaviest states (3 ++ and above) have slightly shorter lifetimes due to having more available decay channels with larger mass splittings for the off-shell Higgs.
One of the most important characteristics of the LHC signatures of dark glue showers is the distribution of glueball lifetimes.Depending on the fundamental parameters of the model, the predictions range from semivisible jets (all the decays that are visible to ATLAS/CMS are prompt) to emerging jets (the glueball decays occur within the ATLAS/CMS detector) to the long lifetime regime (all glueballs escape the main detectors).For small m 0 and large M/y, all species are sufficiently longlived that a dedicated long-lived particle experiment such as MATHUSLA [84][85][86] can significantly extend sensitivity beyond main detector searches due to its large volume.
There are also regions where one kind of jet signature dominates over another, or a mixture of both strategies is potentially viable.Given that a generic strongly-coupled sector can have a spectrum of many hadrons with a broad hierarchy of lifetimes (as in the SM itself), an optimal search could incorporate methods from the semivisible jet, emerging jet, and external LLP detector strategies simultaneously.
In the remainder of this section, we discuss two production mechanisms for glueballs at the LHC: through the Higgs and through a new heavy Z ′ .For each mechanism, we show the parameter space relevant for semivisible or emerging jet searches, and we discuss predictions for two different classes of glueball collider phenomenology.For Higgs portal production, we make predictions for glueball decays that could be observed in the proposed MATHUSLA experiment, considering glueball production and decays through a Higgs portal within FSUSY and FTH as discussed in Sec.3.2.This will supersede the rudimentary MATHUSLA sensitivity estimates for neutral naturalness presented in [87].
For the Z ′ production, we show that this model could yield a good benchmark for semivisible jet and emerging jet searches.One phenomenological parameter used in the studies of semivisible jets is r inv , the average ratio of the number of dark hadrons that are stable on collider scales compared to the total number of total dark hadrons produced.In current searches, one takes r inv as a simplified model-like input parameter and models the distribution of the invisible fraction of hadrons as Poissonian.Our ability to model hadronization and decay of dark glueballs allows us map parameters in the fundamental description onto a prediction for r inv event by event.Thus, we provide a theoretically motivated range of r inv distributions to consider for future semivisible jet searches.

Dark Glueballs via Higgs Production
In this section, we discuss the signatures of dark glueball showers produced via the Higgs portal, which we outline in Sec.4.1.1.Fig. 9 shows fractions of dark glueball events that could possibly give rise to an emerging jet signature, as well as fractions of events that could possibly have a semivisible jets signature with no displaced decays.These plots reveal how different regions of parameter space motivate different  For the emerging jet fractions, events were required to have at least one glueball decay within the CMS tracker with transverse displacement of at least 50 mm [21].For the semivisible jet fraction, events were required to have at least one glueball escape the tracker, at least one prompt glueball decay within the tracker, and no glueball decays within the tracker with transverse displacement > 50 mm.combinations of main detector search strategies depending on the glueball lifetime hierarchy.Therefore, the below sensitivity analysis for MATHUSLA will demonstrate where in parameter space the dedicated LLP strategy has reach beyond the main detectors and where these strategies have potential overlap.

Higgs Production
To model Higgs production of dark glueballs we simulate gluon-gluon fusion and VBF in MadGraph5 amc@nlo [108] + Pythia 8 [83].Gluon fusion is implemented via the effective ggh operator, with jet matching for up to one extra hard jet and slight event reweighting to reproduce the NLO+NNLL Higgs p T spectrum computed by HqT 2.0 [109,110].As discussed in [44], the Higgs-to-dark gluon branching ratio can be found by a rescaling of the SM Higgs-to-gluon branching ratio of 8.5% [111,112].

Dark Glueballs at MATHUSLA
In Fig. 10, we show sensitivity curves for decays within the 100 m × 100 m × 25 m MATHUSLA decay volume as specified in [86], assuming an integrated luminosity of 3 ab −1 at √ s = 14 TeV.The contours we show correspond to 4 decays in MATH-USLA's decay volume, illustrating the exclusion reach in the absence of backgrounds, which is expected for LLP decays to high multiplicities of SM hadrons.The experimental bound on the Higgs-to-invisible branching ratio of 18% [113] excludes the parameter space of top partner masses below the range shown in the plots.We account for uncertainty in the various heavy glueball lifetimes by varying the cor-  We take 4 events within the decay volume as the exclusion limit.The top plots show exclusive decays of the two lightest glueball species, and the bottom is inclusive of all species.The dashed contours reflect uncertainties due to variation of both the hadronization benchmark and the decay matrix elements.The inclusive plot also shows the previous estimate from [87] based on the simplifying conservative assumption of two-body exotic Higgs decays h → 0 ++ 0 ++ only.responding decay constants independently by a factor of 2 in each direction.Our uncertainty bands also include the variation obtained by running simulations with the different hadronization benchmarks introduced above in Table 1.
A striking feature of these results is the importance of including the heavier glueball species, and the resulting dramatic increase to MATHUSLA's estimated sensitivity in neutral naturalness parameter space.Since the heavier glueballs have longer lifetimes than the 0 ++ , MATHUSLA is able to probe an entirely complementary mass regime with heavier glueball decays, extending its reach up to m 0 ∼ 50 GeV compared to the ∼ 20 GeV maximum probed by the 0 ++ alone.
Previous studies of dark glueball phenomenology [44,87] made the conservative simplifying assumption that only two glueballs were produced in exotic Higgs decays, some fraction of which was h → 0 ++ 0 ++ , the only channel assumed to be observable.This was necessitated by the absence of a realistic simulation framework for glueball production.Our work significantly improves on these previous sensitiv-ity estimates by including the dark gluon shower and Lund string-inspired glueball hadronization, allowing for both the higher glueball multiplicity in each Higgs decay and the contribution of heavier, more long-lived glueball states to be systematically taken into account.This leads to the improved projections for the inclusive total reach of MATHUSLA to all glueball decays, shown in the bottom panel of Fig. 10.It is interesting to note that this updated inclusive MATHUSLA reach therefore not only includes the long 0 ++ lifetime regime below 20 GeV, but also exceeds, or is at least comparable to, the total projected coverage of main detector searches relying on LLP decays in the tracker or muon system for m 0 ≲ 50 GeV, computed with the above simplifying two-body-decay assumption [44].While the main detector search sensitivities would be expanded due to increased glueball multiplicity in our updated simulations, the inclusion of heavier glueballs would have a much smaller effect than it did for MATHUSLA, since the main detector is most sensitive to short lifetimes.While the detailed study of main detector sensitivities to dark glueballs is an important subject of future study with our updated simulation framework, this nonetheless already suggests that MATHUSLA's LLP sensitivity may dramatically enhance new physics coverage in a large region of dark glueball parameter space.
Our results also motivate further study into the properties of the heavier species.In particular, lattice computations to determine the decay matrix elements would reduce uncertainty in the glueball lifetimes.In the regions of parameter space where the heavier species dominate decays in MATHUSLA, the uncertainty due to lifetime variation is larger than that due to the hadronization benchmark variation.
Note that the final states of 2 ++ decay always include a 0 ++ , and the region of parameter space where the 2 ++ dominates decays in MATHUSLA is also where the 0 ++ has short O(cm) lifetimes, see Fig. 8. Therefore, given the cm-scale tracking resolution of the MATHUSLA experiment, the 2 ++ decay can be treated as a single vertex.This region of parameter space is also interesting because any 0 ++ produced would decay within CMS.These can be searched for with dedicated searches using CMS detector information alone [44,114] (though with significant signal penalty due to trigger limitations) or a combined MATHUSLA-CMS search if MATHUSLA provides a trigger signal to CMS [115].In the latter case, simultaneous reconstruction of the 0 ++ and 2 ++ decay would allow a detailed characterization of the dark sector and provide strong evidence that the newly discovered LLP states are in fact dark glueballs.

Dark Glueballs via Z ′ Production
In some of our parameter space, the lightest dark glueballs can decay promptly while the rest are either stable or very long lived.This would lead to LHC events where the visible jet transverse momentum ⃗ p J T and missing transverse momentum ⃗ p miss T are aligned, which is the characteristic property of so-called semivisible jets [15].To define a benchmark for future semivisible jet searches, we consider the simplified m 0 (GeV)  For the emerging jet fractions, events were required to have at least one glueball decay within the CMS tracker with transverse displacement of at least 50 mm [21].For the semivisible jet fraction, events were required to have at least one glueball escape the tracker, at least one prompt glueball decay within the tracker, and no glueball decays within the tracker with transverse displacement > 50 mm.signal model used in the recent CMS semivisible jet search [19].This search assumed the resonant production of a Z ′ mediator that decayed to dark sector quarks, which showered and formed dark hadrons that decayed to SM quarks.We retain the Z ′ mediator, but we work in the region of parameter space that produces dark glueballs, and we introduce the Higgs portal to facilitate the glueball decays.Further details of the production mechanism via a heavy Z ′ are outlined below in Sec.4.2.1.In Fig. 11, we show the fractions of dark glueball events that could give rise to a emerging jet and/or a semivisible jet signature through Z ′ production.Semivisible jet searches may be able to probe the large m 0 regime, while emerging jet search strategies may be able to probe parameter space that includes lower m 0 values.The actual ATLAS or CMS sensitivity to these search strategies requires detailed modeling of the emerging or semivisible jets including SM backgrounds, which is beyond the scope of this paper.

Z ′ Production
The signal model features a Z ′ that couples to both SM quarks and dark sector quarks Q D charged under the dark QCD with confinement scale Λ D .This allows for dark quark pair production in LHC collisions pp → Z ′ → Q D QD .In the parameter regime analyzed by the CMS search, the dark sector quarks have mass M Q < Λ D ≪ M Z ′ , which hadronize into jets of dark mesons (bound states of quark-anti-quark pairs).For the semivisible jet benchmark introduced here, we instead consider the quirk-like regime [116], with Λ D ≪ M Q ∼ M Z ′ /2.This implies that Q D QD pair production via the Z ′ resonance results in a quirk bound state.The Q D QD pair are connected by an oscillating flux tube, which de-excites by radiating glueballs (and angular momentum) before the dark quark pair annihilates into dark gluons.Since the de-excitation sheds angular momentum, the final annihilation is anticipated to be dominated by the s-wave.We also assume that the dark quarks couple to the SM Higgs with a Yukawa coupling y.Integrating out the dark quarks generates the Higgs portal operator, which we assume provides the dominant channel for the dark glueball decays.
The dynamics of quirk de-excitation via glueball emission are not well understood, but as a naïve first guess, we assume the glueball radiation from de-excitation is highly subdominant compared to the glueballs produced in the ultimate s-channel annihilation of the Q D QD pair.This can be guaranteed by setting M Q just below M Z ′ /2, where a tiny mass difference is required to allow emission of a single glueball to shed the quirk's orbital angular momentum.This model technically contains both Higgs and vector portals to the dark sector.In practice however, the vector portal dominates dark quark pair production for M Z ′ up to a few TeV, while the Higgs portal dominates glueball decays with lifetimes shown in Sec.3.2.The vector portal glueball decays are phase space suppressed, since they induce four body decays compared to the two body decays that are induced by the Higgs portal.Additionally, the vector portal decays have lower rates due to the higher dimension of the corresponding effective operator.
The existence of the Higgs portal accommodates the heavy quarks being vectorlike doublets under SM SU (2) L × U (1) Y , but one can also consider the same effective operator in an FTH-like scenario where the dark quarks are SM singlets that couple to a scalar that mixes with the SM Higgs boson, resulting in an analogous M Q /y that sets the glueball lifetimes in combination with m 0 .Whether the UV model has SU (2) L doublet quarks that do not get all of their mass from SM electroweak symmetry breaking, or FTH-like quarks whose M Q /y is not fixed by neutral naturalness considerations, we will simply vary M Q ≃ M Z ′ /2 and y independently for the purpose of studying the semivisible jet signal at the LHC.

Dark Glueball Semivisible/Emerging Jets
For this study, we assume the model described in Sec.4.2.1.We take a benchmark where the Z ′ has mass M Z ′ = 3 TeV, and the dark quarks have M = M Q ∼ M Z ′ /2.We choose M/y = 4.5 TeV, which fixes the dimension-6 glueball lifetimes for any choice of m 0 .To generate events, we use MadGraph5 amc@nlo [108] with a Z ′ model [117,118] to simulate pp → Z ′ production at a 14 TeV proton collider.We run our dark shower and glueball hadronization algorithm as though the Z ′ were a heavy scalar decaying to two dark gluons, which models s-wave quirk annihilation as , where r inv is the fraction of dark hadrons that are invisible to the semivisible jet reconstruction.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.Means µ and standard deviations σ are displayed, with uncertainties corresponding to hadronization variations.discussed in Sec.4.2.1. 10We take events with Z ′ production and decay to glueballs as the hard process and pass them to Pythia version 8.307 [83] to handle SM QCD initial-and final-state radiation and jet clustering with the FastJet version 3.4.0plugin [119].Following the procedure from the CMS search in [19], we use the antik T jet clustering algorithm [120] with jet radius R = 0.8.SM final states in each event are characterized as "invisible" if they are neutrinos or they have a glueball ancestor that decayed outside of the cylinder spanned by the CMS tracker [121].We consider all other SM final states as "visible" and cluster them into jets.This is a highly simplified picture of the detector and semivisible jet reconstruction, but we performed the same analysis using the central hadronic calorimeter as the border between visible and invisible states and found qualitatively similar results.The invisible SM final states, as well as any stable glueballs, contribute to the missing transverse momentum ⃗ p miss T .Most of the glueball species decay by emitting a lighter glueball, resulting in cascade decays.However, only the primary glueballs (i.e.those produced from hadronization) contribute to r inv .In order to understand the role of displaced ver- , where r dec is the distance of glueball decay vertices within the CMS tracker to the IP.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.tices, we tracked the distances r dec between the interaction point and the decay vertices of glueballs that decayed within the tracker.We also computed a few useful observables for semivisible jet searches and describe them further App.C. The distributions of these other observables have the expected qualitative form, so here we focus on the novel results of r inv and r dec shown in Fig. 12 and Fig. 13, respectively.
We see that the average value of r inv ranges from ∼ 0.45 to ∼ 0.82.This spread of average r inv demonstrates that a constraint on r inv from a semivisible jet search could potentially be recast as a constraint on our model's microscopic parameters.It would also be interesting to investigate the extent to which the shape of the r inv distribution predicted here would impact limits set by existing analyses.As m 0 increases, the mean r inv approaches ∼ 0.45, and we expect from Fig. 5 that about half of the glueballs produced are 0 ++ .Therefore, this behavior of r inv shows that as m 0 increases through this range, all the 0 ++ glueballs decay within the tracker while all heavier species tend to escape.The r dec plots emphasize the importance of displaced vertices in this analysis.There are two ways for a dark sector to generate missing momentum aligned with a jet: the jet contains states with long lifetimes compared to the detector scale as well as states that decay promptly, or the jet contains states with lifetimes comparable to the detector scale, allowing a portion of them to decay in the detector and leave displaced vertices.The former case is the prototypical semivisible jets scenario, while the latter is closer to an emerging jets signal.With our dimension-6 glueball decays through a Higgs portal, these two cases overlap.Depending on the region in parameter space, the 0 ++ may decay promptly or have cm to m lifetimes while the heavier states are relatively long-lived.Alternatively, the 0 ++ may decay promptly while only a subset of the heavier species leave displaced vertices.We leave the interesting task of developing an optimal search strategy that takes advantage of this class of signals to future work.

Conclusions
Confining dark sectors appear as a component of a large class of possible BSM scenarios with a wide range of possible signatures and very broad theoretical motivation.The case of a pure Yang-Mills dark sector, corresponding to N f = 0 QCD, is an important representative of this class, but one whose study has been hampered by our ignorance of pure-glue hadronization into dark glueballs.While first steps to constrain the possible range of hadronization outcomes based on the perturbative gluon shower, energy conservation, and the large m 0 /Λ D mass gap were taken in [81], the absence of a realistic hadronization model left large theoretical uncertainties, especially for exclusive production rates of individual glueball species, which determine most collider observables.This motivates developing a more sophisticated phenomenological parameterization of the non-perturbative physics, which allows us to make more quantitative predictions for the final state from dark glue sector showers.
In this work, we present the first implementation of a color string dynamics inspired hadronization model for dark glueball fragmentation.We borrow from the Lund string model to parameterize how color-singlet flux tubes produced by the dark gluon shower fragment into dark glueballs.The specific algorithm is quite simple, and local in the sense that it assumes each combination of fragmenting string pieces chooses democratically from all kinematically available glueball species.This makes one of our main results all the more remarkable: the relative multiplicities of different glueball species approximately follows the theoretically expected thermal distribution with T had ≃ T c ≃ Λ D , independent of the chosen fragmentation function and with a weak dependence on the shower cutoff scale.Also, our algorithm's intuitive geometric picture of self-intersecting color flux rings suggests our approach may capture the relevant dynamics of glueball fragmentation.
Our shower and hadronization model has a handful of parameters beyond the physical glueball mass and initial center of mass energy: the shower cutoff scale, the choice of fragmentation function, and that function's two parameters.We suggest three sets of benchmarks that span a broad space of our model's physically reasonable predictions for different possible values of these nuisance parameters, see Table 1.
We apply our new simulation framework to explore the potential reach of emerging and semivisble jet searches, finding both search strategies would probe distinct but overlapping regions of dark glueball parameter space.We then focus on two important phenomenological demonstrations: the study of glueball decays in the proposed MATHUSLA detector and a realistic benchmark model for semivisible/emerging jets.
For dark glueballs decaying in MATHUSLA, we focus on the simplified glueball parameter space motivated by theories of neutral naturalness like the Fraternal Twin Higgs or Folded Supersymmetry, where the dark QCD sector is coupled to the SM via the Higgs portal.We find that including the multiplicity-enhancing effects of the shower and the realistic full spectrum of glueball species produced from our hadronization model dramatically increase the projected sensitivity of MATHUSLA compared to earlier simplified estimates, see Fig. 10.
We also considered a simplified model of a Z ′ coupled to a dark QCD in the quirklike regime that lead to the production of dark glueballs via resonant Z ′ production.This model has broad and potentially overlapping regions of parameter space that can possibly have emerging jet and semivisible jet signatures.The main result we obtained using our simulation of the pure-glue shower and hadronization in the Z ′ model was finding the distributions of r inv that arise due to the pure-glue theory's multiplicity of glueball states with potentially widely separated lifetimes, see Fig. 12.This demonstrates how a semivisible jet search can yield realistic constraints for pure-glue dark sectors.
We hope that some version of our approach can be incorporated into the Pythia Hidden Valley Module.For future studies, we suggest further lattice calculations of glueball decay matrix elements, which will reduce systematic uncertainties in glueball lifetimes relevant to the long-lived particle regime.Our hadronization model could also be improved or generalized in various ways, such as accounting for parity and charge selection rules, and generalizing the number of dark colors beyond 3 (though fully characterizing the glueball mass spectrum for N c ̸ = 3 would again require additional lattice studies).A more sophisticated implementation might explicitly simulate a closed oscillating classical string in position space.This would improve the IR/collinear unsafety inherent in our approach, which relies on discretizing rings into string pieces in a way that explicitly depends on the number of gluon splittings during the shower.
A reliable glueball Monte Carlo also enables new cosmological studies of these dark sectors.Previous analyses derived constraints from avoiding overproduction of surviving stable glueball states [64,65,73,75,122,123] or late-time decays modifying big-bang nucleosynthesis and the cosmic microwave background [74].These analyses assumed that the glueball relic densities originated in the dark QCD phase transition, which itself is not well understood.However, it is possible for glueball densities to receive important contributions from late-time decays, annihilations, or other entropy injections.For these scenarios, our shower and hadronization model could supply new predictions.The same can be said for models where glueball production plays an important role in the production of cosmic rays [71,124].
The most immediate application of our work will be to enable many new detailed collider studies of N f = 0 QCD dark sectors.For example, including factors such as detector effects is required in order to quantitatively establish the distinction between the emerging and semivisible jet regimes, as well as the realistic sensitivity of the ATLAS and CMS main detectors to dark glueball decays.This work also provides us with a tool we can use to further develop search strategies that are insensitive to dark hadronization uncertainties [93], perhaps relying on some aspects of jet substructure and that could even incorporate machine learning, see e.g.[18,28,[125][126][127][128][129][130][131][132].In general, understanding the detailed phenomenology of these dark sectors will help us design the searches that could lead to the next discovery beyond the Standard Model.
λ is the seed for glueball emission. 11The two string pieces connected by this vertex are (at least) the pieces that will be converted or "fragmented" into a glueball.If the two pieces have an invariant mass less than m 0 , then one of their nearest neighbors (selected arbitrarily) will also be added to the list of fragmenting pieces.If this is still not enough invariant mass, the nearest neighbor in the other direction (still in color-connection space) will fragment as well, and so on until we have a selection of string pieces with invariant mass at least m 0 .
Next, we select the glueball's species.The invariant mass of the fragmenting string pieces is an upper bound on the mass m G of the selected species.We perform a weighted random selection from the kinematically available species, where the weight 2J + 1 accounts for the higher multiplicity of a state with spin J.A similar approach is used for SM meson species selection in Pythia.For example, a u d mode could be either a pion or a ρ meson.The choice is determined in Pythia by incorporating a spin-dependent weight into the random selection.As a correction to better fit the data, Pythia additionally includes suppression of heavier species beyond the naïve 3:1 expectation for the vector-to-scalar weight. 12As an alternative approach, one could plausibly select the species before the fragmenting string pieces and gather string pieces with invariant mass at least m G .However, that approach causes the species distribution to acquire a strong dependence on the ratio of the invariant mass of the shower, which is distinctly different behavior from what we expect in QCD as discussed in Sec.2.3.After determining the glueball species, we must specify its momentum.First, we choose the direction of the three momentum pG to be along the momentum of the fragmenting string pieces in the ring's rest frame.This decision is admittedly not Lorentz covariant, but this pG is the best indication of a preferred direction of the fragmenting system.If all of the string pieces in the ring are fragmenting, then pG is selected randomly and isotropically in the ring's rest frame.Then, again inspired by the Lund string model, the glueball takes a fraction z of the fragmenting pieces' light-cone momentum p ±pieces = E pieces ± |⃗ p pieces |, so that the glueball's light-cone momentum is which is a relation that is invariant under boosts along pG .With the glueball's direction and light-cone momentum fixed, and imposing the on-shell condition, the

B Benchmark Parameter Variations
Here, we show variations in the species production rates and glueball momentum distributions resulting from our different benchmark parameters.The momentum distributions in Fig. 14 provide an intuitive demonstration of what we mean by glueball "hardness."As discussed in Sec.2.4, there is a straightforward relation between our hadronization algorithm's nuisance parameters and the glueballs' tendency to be produced with smaller or larger momentum.In this sense, the soft and hard variations of our suggested benchmarks are meant to provide the extremes of our algorithm's possible sensible outputs.
As seen in Fig. 15, as the shower cutoff scale (parametrized by c) increases, the best-fit T had also increases.This behavior is as expected because a higher shower cutoff scale means the dark gluons have fewer opportunities to branch, leading to higher-mass string pieces during fragmentation and weaker suppression of heavy species production.As in Figs. 4 and 1, measured in the rest frame of the dark gluon shower.Exclusive distributions of the two lightest species are shown, as well as the inclusive distribution.As expected, glueballs from "harder" parameter variations tend to have larger momentum.

C Additional Semivisible Jet Distributions
In addition to r inv and r dec , we computed three observables that [19]  and the two highest-p T jets.The m T distribution is essentially cutoff by the mass of the Z ′ , and R T and ∆ϕ min can be used to help cut out background.We found that the ∆ϕ min and R T distributions were fairly consistent across the model parameters we simulated, so we show a few representative examples in Fig. 16.The m T distribution changed more noticeably, so we show more model parameter variations in Fig.
Partons at the end of the shower with N c → ∞ have unique color and anticolor labels.String pieces are randomly reassigned a new color, now restricted to N c = 3 choices.
between color/anti-color pairs are swapped if this reduces the string-length λ.

Figure 1 :
Figure 1: A sketch of QCD-based color reconnection.Color and anti-color labels are represented by filled and empty (semi-)circles, respectively.Quarks q and antiquarks q form the ends of strings, while gluons g are kinks in the string.String pieces span between the color/anti-color pairs of individual partons and are labeled with displayed colors (colors have numerical labels so that gray-scale versions contain the same information).Each parton is displayed with a fixed location in an abstract color-connection space where the displayed lengths of the string pieces correspond to the string-length λ. Between 1b and 1c, the pairs of blue (1) and green (3) connections are swapped.
Gluons at the end of the N c → ∞ shower.
Reassignment of string piece colors for N c = 3.

Figure 2 :
Figure 2: A sketch of QCD-based color reconnection for a system of only gluons, as in Fig. 1. Between 2b and 2c, the pair of blue (1) connections are swapped.
(a) The vertex (green) joining the string pieces with the largest string-length begins the fragmentation.(b) A minimal set of string pieces (blue) with total mass ≥ m 0 nearest to the vertex is selected to fragment into a glueball.(c) The glueball (red circle) takes a fraction of the fragmenting pieces' momentum.The remaining momentum is distributed among two new string pieces (blue).

Figure 3 :
Figure 3: A visual depiction of our glueball hadronization mechanism.Color-singlet rings are shown as polygons whose edges are Lund string pieces.The criterion for selecting the seed vertex (green) is chosen to allow for the most rapid reduction in overall string length with each step.Glueball emission (depicted by the red circle in 3c) is conceptualized as pinching off from the fragmenting string pieces (blue in 3b).

Figure 4 :
Figure4: Thermal model fit quality (top) and corresponding best fit hadronization temperature in units of Λ D (bottom) using the beta distribution (left) and the LSFF (right) for points in the plane of the mean µ z and standard deviation σ z of the 0 ++ fragmentation function.The shower center of mass energy is 125 GeV, the lightest glueball mass is m 0 = 10 GeV, and the shower cutoff parameter is set to the default c = 1.8.The black contours indicate an upper bound on σ z for a given µ z .For the beta distribution, this is the least upper bound for parameters where f β (0) and f β (1) are finite.For the LSFF, σ z does not saturate the bound.

Figure 7 :
Figure 7: Average glueball multiplicity as a function of E CM .The points show the output of our algorithm, and the solid lines show the QCD prediction in Eq. (2.7) normalized to match the Monte Carlo at the largest E CM .

Figure 8 :
Figure 8: Contours show log 10 (cτ /m), where cτ is the mean decay length of the glueball in the space of the lightest glueball mass m 0 and top partner masses assuming folded supersymmetry m tFSUSY or fraternal twin Higgs m tFTH .The top plots show the two lightest species, and the bottom plots show representative examples of heavier species.

Figure 9 :
Figure9: Fractions of dark glueball events for the Higgs production scenario satisfying necessary but not sufficient conditions to produce emerging jet signals (left) or semivisible but not emerging jet signals (right).For the emerging jet fractions, events were required to have at least one glueball decay within the CMS tracker with transverse displacement of at least 50 mm[21].For the semivisible jet fraction, events were required to have at least one glueball escape the tracker, at least one prompt glueball decay within the tracker, and no glueball decays within the tracker with transverse displacement > 50 mm.

Figure 10 :
Figure 10: Sensitivity curves for glueball decays in MATHUSLA in the space of the lightest glueball mass m 0 and top partner masses in the fraternal twin Higgs model m tFTH or folded supersymmetry model m tFSUSY .We take 4 events within the decay volume as the exclusion limit.The top plots show exclusive decays of the two lightest glueball species, and the bottom is inclusive of all species.The dashed contours reflect uncertainties due to variation of both the hadronization benchmark and the decay matrix elements.The inclusive plot also shows the previous estimate from[87] based on the simplifying conservative assumption of two-body exotic Higgs decays h → 0 ++ 0 ++ only.

Figure 11 :
Figure11: Fractions of dark glueball events in the Z ′ production scenario with m Z ′ = 3 TeV satisfying necessary but not sufficient conditions to produce emerging jet signals (left) or semivisible but not emerging jet signals (right).For the emerging jet fractions, events were required to have at least one glueball decay within the CMS tracker with transverse displacement of at least 50 mm[21].For the semivisible jet fraction, events were required to have at least one glueball escape the tracker, at least one prompt glueball decay within the tracker, and no glueball decays within the tracker with transverse displacement > 50 mm.

Figure 12 :
Figure 12: Distributions of r inv for various values of the lightest glueball mass m 0 in the Z ′ production model with m Z ′ = 3 TeV andM Q ∼ M Z ′ /2, where r inv is the fraction of dark hadrons that are invisible to the semivisible jet reconstruction.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.Means µ and standard deviations σ are displayed, with uncertainties corresponding to hadronization variations.

Figure 13 :
Figure 13: Distributions of r dec for various values of the lightest glueball mass m 0 in the Z ′ production model with m Z ′ = 3 TeV andM Q ∼ M Z ′ /2, where r dec is the distance of glueball decay vertices within the CMS tracker to the IP.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.
5, Figs. 14 and 15 were generated with m 0 = 10 GeV and dark shower center of mass energy of 125 GeV.

Figure 14 :
Figure 14: Distributions of |⃗ p |/m 0 for the three sets of benchmark parameters listed in Table1, measured in the rest frame of the dark gluon shower.Exclusive distributions of the two lightest species are shown, as well as the inclusive distribution.As expected, glueballs from "harder" parameter variations tend to have larger momentum.

Figure 16 :
Figure 16: Distributions of ∆ϕ min and R T for various values of the lightest glueball mass m 0 for the semivisible jets scenario.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.

2 T = m 2 JJ + 2|⃗ p miss T | m 2
used in their search.The transverse mass m T is given bym JJ + |⃗ p T,JJ | 2 − |⃗ p T,JJ | cos(ϕ miss JJ ) , (C.1)where the two highest-p T jets have total momentum p JJ with corresponding invariant mass m JJ , and ϕ miss JJ is the azimuthal angle between ⃗ p T,JJ and ⃗ p miss T .The other two observables are R T = |⃗ p miss T |/m T and the minimum azimuthal angle ∆ϕ min between ⃗ p miss T

17
. Since R T depends on m T but has significantly weaker dependence on the model, there must be a compensating change in |⃗ p miss T | as m T changes.All of these plots were generated with M/y = 4.5 TeV.

Figure 17 :
Figure17: Distributions of m T for various values of the lightest glueball mass m 0 for the semivisible jets scenario.Solid histograms come from using the default hadronization benchmark, and dashed histograms come from the soft and hard variations.
Relative production rates of primary SM hadron species, normalized to number of spin degrees of freedom, following decay of a 1 TeV scalar to two gluons in Pythia with electroweak interactions turned off.The grouping is determined by identifying the heaviest valence quark flavor within the hadron, and each group is fit to its own Boltzmann distribution (solid lines) by maximizing the coefficient of determination R 2 .The contribution to the fit from hadron/anti-hadron species pairs are averaged to avoid fitting the same masses twice.The Monte Carlo statistics are sufficiently high that the excesses of heavy hadrons compared to the fits are not due to random fluctuations.

Table 1 :
D (p T min ) µ z σ z T had /Λ D Suggested benchmarks to set the hadronization scale and fragmentation function nuisance parameters.Also shown are the corresponding values of the dark sector coupling α D evaluated at the shower cutoff scale p T min = c Λ D , the mean µ z and standard deviation σ z of the 0 ++ fragmentation function, and the best-fit T had in Eq. (2.5) for the relative species production rates.

Table 2 :
Table of masses and decay channels for each glueball; h * indicates an offshell Higgs.