Spheres to jets tuning event shapes with 5d simplified models

Hidden sectors could give rise to a wide variety of events at the LHC. Confining hidden sectors are known to engender events with a small number of jets when they are weakly-coupled at high energies, and quasi-spherical soft unclustered energy patterns (SUEPs) when they are very strongly-coupled (large ‘t Hooft coupling) at high energies. The intermediate regime is murky, and could give rise to signals hiding from existing search strategies. While the intermediate coupling regime is not calculable, it is possible to pursue a phenomenological approach in which one creates signals that are intermediate between spherical and jetty. We propose a strategy for generating events of this type using simplified models in extra dimensions. The degree to which the event looks spherical is related to the number of decays produced near kinematic threshold. We provide an analytic understanding of how this is determined by parameters of the model. To quantify the shape of events produced with this model, we use a recently proposed observable — event isotropy — which is a better probe of the spherical regime than earlier event shape observables.


Introduction
Data taken at the Large Hadron Collider (LHC) has been essential to confirm Standard Model (SM) predictions and perform precision measurements [1,2]. Extensive efforts have been made to search for phenomena beyond the Standard Model (BSM) as well, such as supersymmetry, dark matter, extra dimensions, and many others. Thus far, no BSM physics has been observed. The lack of evidence for new physics suggests that the simplest scenarios are not good models of nature.

JHEP05(2021)096
However, one should worry that more complicated signals might have evaded current LHC search strategies. A large class of unusual BSM signals arise within the "hiddenvalley" (HV) scenario, in which the SM is weakly coupled to a hidden sector of selfinteracting particles neutral under the SM gauge groups, at least one of which decays visibly [3][4][5][6]. Similarly unusual signals arise in theories of "quirks" [3,[7][8][9], and also in the context of unparticles [10][11][12] which in the presence of a mass gap may produce hiddenvalley phenomenology [13]. Since hidden sectors can be arbitrarily complex, many models of this type are poorly constrained by LHC searches, and since indirect limits on hidden sectors are often very weak, there are few other constraints. But it is impossible to search through the full space of HV models, or of models in other general scenarios. Instead, the practical approach to discovering a new signal is to carry out searches for distinct and parametrizable signatures, and to make this possible, models that produce these signatures must be developed. Among the unusual signatures identified in HV models to date are soft, unclustered energy patterns (SUEPs) [13][14][15]; lepton jets [16,17]; emerging jets [3,18]; semi-visible jets [3,[19][20][21]; and dark jets with unusual substructure [22,102].
General, flexible searches for new physics at the LHC require us not only to widen the range of models and signatures, but also to expand the tools available for data analysis, both offline and at the trigger stage (where fewer than 1% of LHC collisions are recorded). Search strategies for general hidden sectors are confounded by the many free parameters, including the masses, couplings, and lifetimes of new particles. The overall event shape of different scenarios can also take many forms. This motivates the development of diagnostic tools to characterize anomalous event shapes which are unlikely to arise from Standard Model processes.
In this paper we will introduce a class of simplified models that produce a range of new signatures, and our goal will be to characterize them using event shape observables. It would be premature to consider how to search for these signals at the LHC, as we should first understand the signatures themselves. For this reason we focus our attention on an idealized situation: a pure signal at an e + e − collider of the future with no background.
There exist several well-known observables that characterize the shapes of events at e + e − colliders. A commonly used infrared and collinear (IRC) safe event shape observable is thrust, defined as [23][24][25] It has a range T ∈ [0.5, 1], where T = 1 corresponds to two back-to-back particles, and T = 0.5 is an isotropic radiation pattern. While thrust has provided essential insight on the perturbative nature of QCD, it is most sensitive to event shape deviations from the two-particle dijet configuration and has less sensitivity in the quasi-spherical regime (see figure 8 of [26]). To complement such standard observables, we also make use of a recently proposed event shape observable, the event isotropy [26]. Event isotropy is defined using the energy mover's distance (EMD) [27,28]. Given two radiation patterns of massless particles P , Q, the EMD is the minimum work necessary to rearrange P into Q. A radiation pattern is defined as a set of particles, each of which are specified by their position on the unit sphere and fraction of the total energy. To reorganize values than the corresponding plot in ref. [26].
the pattern P into Q, we construct a transport map f ij that tracks the total energy fraction moved from position i to position j. The total work done in the rearrangement can be written as a sum over the distance d ij from i to j weighed by the fraction of energy moved f ij . The EMD is the minimum work for all possible rearrangements of P to Q: (1. 2) The distance measure we use in this paper is forn i the unit vector proportional to the three-momentum of the element p i , etc. Note this differs from [26], where the distance measure was proportional to d ij ∼ 1 − cos θ ij . See appendix B of that paper for a discussion of the different distance measures.
The ideal event isotropy of an event would be its EMD to a uniform radiation pattern U of equal total energy: (1.4) which we would take as a uniform spherical distribution for e + e − colliders. However to make computation times practical we calculate the EMD to a uniformly tiled, high multiplicity sphere generated with HEALPix [29]. Specifically, we will use multiplicity 192 for the tiling, and we denote the event isotropy variable by I sph 192 . Note that event isotropy is defined as a distance to an isotropic event. Thus, highly isotropic events have low values of I sph , despite the name. In particular, I sph is 0 for a spherical event and 1 for a pencil dijet event with two back-to-back momenta. More generally, for an event with k isotropically distributed particles of equal energy, a sphere . (1.5) The minimum value of I sph 192 ≈ 2 192 ≈ 0.1 represents a slight loss of dynamic range, a price one pays for much faster computation.
To illustrate the behavior of this variable, we used Pythia 8.243 [30] to generate, at a center-of-mass energy of 350 GeV, the processes e + e − → qq and e + e − → tt. The resulting event isotropy distributions are shown in figure 1; the values of I sph 192 are higher than the corresponding plot in ref. [26] because of our choice of distance measure eq. (1.3). QCD radiation and radiative-return to the Z boson both reduce the qq isotropy from 1 to approximately 0.8. Meanwhile, since the top quarks are produced near threshold, their six jets are distributed quasi-isotropically, and the distribution peaks near the value 1/3 ∼ 0.58 that eq. (1.5) would suggest. HV extensions to the SM often take the form of confining hidden sectors. However, confinement is compatible with a wide range of event shapes. The 't Hooft coupling λ = α s N c plays a major role in determining the shape of events. In a QCD-like, asymptotically free theory, λ 1 near the confinement scale but runs to be 1 at energies well above Λ QCD . Gluon radiation in this regime is characterized by perturbative showering, in which a hard quark or gluon is dressed with moderate amounts of collinear radiation, leading to a classic QCD jet. However, a near-conformal field theory may maintain λ 1 over a wide range of energies above the confinement scale. In this regime collinear radiation is extremely rapid, and all hard partons lose their energy; only soft physics survives. Such a theory will produce events which, because of greatly enhanced radiation [31], are spherically symmetric in the extreme λ → ∞ limit [13,32,33].
Although one can do reliable computations for λ 1 (where field-theory perturbation theory is valid) and sometimes for λ 1 (where gauge/string duality furnishes us with an alternative perturbation expansion), BSM physics could fall in the intermediate regime. This regime is poorly understood, as there are no methods for detailed calculation at intermediate λ.
Because of this obstacle, we pursue a more pragmatic approach. In this paper we seek to develop a procedure for generating events that is physically reasonable and has parameters that allow it to interpolate between jetty and spherical. We hope that a flexible method for producing events with intermediate values of event shape observables may allow the design of analysis and trigger strategies that are sensitive to a broader class of new physics models, even ones whose signals cannot currently be calculated.
One widely-adopted strategy in the search for new physics at colliders is the use of simplified models. These models can abstract away many details of a theory while preserv- 1 We may imagine spreading each particle's energy over a region of area 4π/k, which requires moving the energy a distance of order 1/ √ k. The normalization is obtained by requiring that for k = 2, the extreme dijet limit, I sph 192 = 1. Any non-uniformity in the distribution generally increases the isotropy when averaged over many orientations. For the distance measure used in [26], there is no square root. See appendix A of [26] for more details. ing key elements of the collider phenomenology, and have been designed for a broad range of signals. However, they are intrinsically 'simple': the models include a small number of new light particles and interactions [34][35][36]. While it is beneficial to have fewer parameters, we would also like to consider more complex theories with many interactions and heavier particles, as in Hidden Valleys. There have been studies of simplified models with relatively high final-state multiplicities [19,[37][38][39], as well as studies of specific theories that give signals made of many soft particles [9,14,15]. However, there is a gap which can be filled by a straightforward approach to formulating simplified models that are flexible enough to span a wide range of event shapes.
In this paper, we show that a simplified model within an extra dimension, with a small number of parameters, allows for the generation of a wide range of event shapes. Specifically, we consider a warped extra dimension, in the form of a slice of a 5d AdS space [40]. This is motivated by the AdS/CFT or "gauge/string" correspondence, which relates gauge theories to string theory [41][42][43]. The correspondence suggests that a set of interacting fields (including gravity) on such a finite-size warped extra dimension can be interpreted as the dual of a (large N ) confining gauge theory. (This relation has been made precise for a few purely 4d gauge theories, as in [44,45].) The infinite tower of massive Kaluza-Klein modes (KK modes) in the extra dimension is equivalent to a tower of hadrons of the confining gauge theory. Rather than choosing some ansatz for the masses and couplings of hadrons in a strongly-coupled gauge theory where we cannot do calculations, we choose a small set of masses and couplings in 5d which determine the entire infinite set of masses and couplings in 4d. This reduces the number of arbitrary choices to make, while still allowing enough flexibility to generate a wide range of collider event shapes. This use of AdS as a simplified dual to a Hidden Valley is in the spirit of previous uses of AdS to model low energy QCD [31,[46][47][48][49][50][51]. The fact that KK-mode cascades in a warped extra dimension can produce approximately spherical events was previously examined in [52]. We build on this by constructing a wider range of models, which lead to a wider range of collider events. We also make use of the new tool of event isotropy to obtain an improved characterization of these events.
We emphasize that what we present here is a physical model (morally dual to a field theory at large λ) which can interpolate phenomenologically between the jetty regime (which arises at small λ as it does in QCD) and the quasi-spherical regime (which appears at very high energy at large λ.) Because multiple, qualitatively different interpolations between the two regimes likely exist, our model may have little to do with what one would observe in a real confining gauge theory whose ultraviolet value of λ is varied from small to large. Nevertheless, our approach widens the space of sensible targets for experimenters, and one may hope any search strategies that it inspires may be sensitive to a variety of models, not just the one proposed here, that sit between the jetty and spherical extremes.
We proceed with a brief introduction of simplified models in extra dimensions in section 2. In section 3, we present results from simulations of cascades generated with our models, for a variety of parameter values and interaction terms. We show that these models can accommodate varied distributions of event isotropy. The features of the model are driven by basic properties of the couplings among various KK modes. In simple scenar- ios with one self-coupled bulk field, near-threshold decays are often preferred, while decays with greater available phase space are suppressed. This leads to low-momentum daughter particles with no preferred boost axis, and thus to nearly isotropic events. The degree of isotropy depends on further details, such as the extent to which there is an approximately conserved KK number, and the number of stable KK modes at the bottom of the spectrum. In scenarios with multiple bulk fields, we find that there are "plateaus" in phase space with relatively high decay rates, far from threshold. These lead to less isotropic events. Finally, cases with boundary-localized couplings can have branching ratios determined mostly by phase space, and lead to much less isotropic events. In each case, we explain how the pattern of branching ratios is reflected in properties of the event: thrust, particle multiplicity, the energy distribution of daughter particles, and the new event isotropy observable. The structures that we find in the patterns of couplings among various modes are determined by overlap integrals involving products of three Bessel functions. In section 4, we give an analytic understanding of these integrals. In particular, we show that the overlap integrals can be separated into two terms, one of which can be computed approximately and one of which can be computed exactly. The latter term often dominates, and allows us to obtain a clear analytic understanding of both the regime in which near-threshold decays are preferred and the regime with plateaus of enhanced decays away from threshold. All of the important qualitative features determining the event spectra can thus be extracted from the analytic results. In section 5, we conclude and summarize both forthcoming work and open questions for the future.
A preliminary version of some of our results was reported in section 7.3 of a recent white paper on long-lived particles at the LHC [53]. This also included a comparison to a parton shower algorithm pushed to strong coupling (work of Marat Freytsis), which may be of interest to some readers.
In a companion paper [54], we will provide a more detailed understanding of the event shape observables, and establish that event isotropy captures features of events that are not easily extracted from traditional variables (thrust, eigenvalues of the sphericity tensor, and jet multiplicities).

Extra dimensional simplified models: a brief introduction
One approach to robust searches for new physics at colliders is the use of simplified models. An extensive summary can be found in [55]. Because these models are characterized by effective Lagrangians with only a few new particles, they are not representative of the rich spectra of new particles and decay chains that can arise within hidden sectors. It is unreasonable to select, by hand, the masses and couplings of large numbers of particles. Various approaches to this problem have been chosen in the literature, which we will briefly discuss below in section 2.3. For our purposes, an efficient approach to generating simplified models with a small number of free parameters is to consider theories with an extra dimension containing a small number of bulk fields and interactions, which then produce many modes and couplings in the four-dimensional reduction. This choice has the advantage that gauge/string duality furnishes us with an interpretation of the extra--6 -JHEP05(2021)096 dimensional simplified model in terms of a toy model of a confining hidden sector at large 't Hooft coupling.
Readers familiar with RS models [40] can skim this section: the main message is that we consider scalars with trilinear bulk interactions as a simplified model for the spectrum and essential interactions in the hidden sector.

Spectrum of masses
We will imagine coupling the SM to a hidden sector which consists of states that propagate in at least five dimensions. Such a hidden sector might in principle be dual to a gauge theory via the gauge/string correspondence, and we will often use the language of this correspondence in describing it.
Specifically, let us begin by considering a slice of (4+1)d AdS space (RS1). We will denote the extra spatial coordinate as z. This spacetime geometry is specified by the curvature radius R of the 5d geometry, with metric where z UV < z < z IR . In this paper we take z UV = 0, in order to focus purely on modeling the hidden sector; coupling the sector to the SM may require reintroducing z UV depending on the nature of the interaction between the two sectors. 2 A theory of fields propagating on AdS 5 for z < z IR is often called the "hard-wall" model and has been extensively studied as a model for QCD [31,47,48,50,51]. The dimensionful parameter z IR plays the role of the confinement length scale in pure Yang-Mills theory. Indeed, this type of model is a cartoon of sorts, representing more realistic string constructions that are dual to quasi-conformal field theories which vaguely resemble QCD. More precisely, these field theories are asymptotically conformal at high energy (corresponding to small z) with a continuous coupling constant, and their conformal invariance is broken at a scale Λ that corresponds to z ∼ z IR . In some cases the breaking of conformal invariance is due to confinement. Simple versions of the hard-wall model, like pure Yang-Mills theory, have a mass gap and towers of states, the details depending on the 5d fields that the model contains. We will consider theories of this type below.
For simplicity only, we will consider interacting scalars propagating in the bulk. These could be a subset of fields in a more realistic theory, or could serve as warm-ups for gauge and/or gravity fields. The scalars will satisfy a 5d Klein-Gordon equation with mass M , and can be written as an infinite sum of scalar modes that propagate in the 4d bulk modified by wavefunctions in the fifth dimension: (2.2) 2 SM fields do not propagate in the bulk, because they are not composite states of the hidden sector.
One possibility would be to couple to them to the bulk fields by UV-brane localized interactions, but we will not pursue the details here.  Figure 2. The mass spectra of KK towers for ν = 0, 1, 5, and 10 respectively, in units of 1/z IR starting at the lowest mode. We assume Dirichlet boundary conditions on the IR brane. We highlight two important trends. First, the lowest mass in the tower increases as ν increases. Second, the mass splittings sufficiently high in each tower are approximately equal and independent of ν.
The 5d wavefunctions have the form of Bessel functions: (2. 3) The tower of massive The range 1 ≤ d O < 2 is allowed by unitarity but requires an alternative boundary condition at z UV [57], and will not be considered in this paper. One can estimate the masses of the Kaluza-Klein modes by using the asymptotic expansion of the Bessel function for large (positive real) arguments, In particular, the n th KK mode mass is approximately given by where C is an O(1) number that depends on the choice of boundary condition on the IR brane. For Dirichlet boundary conditions (ψ n | z IR = 0), we have C = −1/2, while for Neumann boundary conditions (∂ z ψ n | z IR = 0), we have C = −3/2. The details of the mass spectrum impose constraints on particle decays; examples of spectra for different ν are given in figure 2. For the Dirichlet case, C = −1/2 implies -8 -

JHEP05(2021)096
that for ν < 1/2, a decay of KK mode n 1 to KK modes n 2 and n 3 of the same field, with n 1 = n 2 + n 3 , is always kinematically allowed. In particular, barring some additional constraints, the only stable mode is n = 1. Conversely, for ν > 1/2, some of these decays are always disallowed, and (to a very good approximation) 1+ (ν + 3 2 )/2 modes are stable against decay. In particular, the mass spectrum for ν = 1 2 is exactly that of a 5d-massless field in a flat extra dimension, m n = πn; as ν → 1 2 from below, the phase space for the decay 2 → 1 + 1 closes off and the second KK mode becomes kinematically stable. Note also that for all such spectra, a decay n 1 → n 2 + n 3 is always forbidden for modes of the same field if n 1 < n 2 + n 3 .
The wavefunction of the nth mode is is determined by requiring that the 4d field φ n (x) be canonically normalized at tree level. The Dirichlet case has a simple normalization: For the Neumann case the closed form is more complicated, but the final approximate expression in eq. (2.9) remains true at large n.
In the case of Neumann boundary conditions, the value of the wavefunction on the IR boundary can be approximated by: This implies that a φ 3 interaction localized on the IR boundary will give rise to couplings of approximately equal magnitude between any three KK modes, a fact that we will make use of in section 3.4.

Interaction terms
Each scalar field in 5d provides a tower of massive particles. To induce the decays among these particles that will create a range of signals, we next turn on interactions among the scalars. In many models, the dominant decays are all two-body, and in such cases, the only important interactions are cubic. We therefore consider cubic couplings of scalar fields in the 5d bulk, where Φ 1,2,3 are potentially different fields with corresponding bulk mass parameters ν 1,2,3 , and c is a coupling constant. 3 We implicitly assume that the unboundedness of this Lagrangian is corrected by higher-order terms which make the theory well-behaved but do not affect decays.

JHEP05(2021)096
In our studies below, we will focus only on the "single field case", where all three fields are the same, and the "two field case" where Φ 2 = Φ 3 . The single-field case captures some features of self-interacting bulk fields such as a dilaton, non-abelian vectors or scalars, or the gravitational field. The two-field case captures features of situations in which a scalar, gauge field, or gravity couples to a second field that carries charge under either a Z 2 symmetry, or perhaps CP, forbidding modes of Φ 2 from being created singly. It also is similar to cases in which the second field is complex and carries a U(1) charge, since the spectrum and decay modes of Φ 2 and Φ * 2 are the same in such a case. The field Φ 1 may have its own cubic interaction, but we will assume here for simplicity that its coupling is relatively small compared to the Φ 1 Φ 2 2 coupling, and so plays an insignificant role in decay chains. The cubic interaction could also be zero, as for an abelian gauge field.
We will denote the i th KK mode of the field Φ n by φ n,i (x). In the 4d effective theory, the 5d interaction translates into an infinite set of couplings among the 4d modes: where the effective couplings of the 4d scalars are determined by the overlaps of the wavefunctions in the extra dimension, (2.13) Here we substituted det g = (R/z) 10 for the metric eq. (2.1) into eq. (2.11).
These 4d coupling constants depend on various dimensionful 5d quantities: c, R, M n , and z IR . However, it turns out that if we write c ijk in terms of the ν n (which depend only on the product M n R) and the dimensionless coupling c 0 ≡ c √ R, then all remaining dependence on R drops out of the equation. Hence, we never need to specify the 5d length scale R to do a calculation. Furthermore, we are primarily interested in branching ratios (rather than total widths), for which the value of c 0 cancels out as well. Consequently, we can set both R and c 0 to 1 for convenience. Meanwhile the scale z IR , corresponding to the confinement scale of a dual field theory, is the only dimensionful quantity that appears in physical measurements. We may express masses and widths of the KK modes, and other dimensionful measurements, in units of z IR , and so we can also set this quantity to 1 if we choose. The only non-trivial parameters left, then, are the 5 dimensional masses, which we express using the dimensionless parameters ν i , which in the context of a gauge/string duality are related to 4d operator dimensions.
As we will see, depending on the choices of ν i and the boundary conditions, the c ijk will often (but not always) respect an approximate KK-number symmetry. Were the extra dimension flat, it would have a conserved KK-number with Neumann boundary conditions (or with periodic boundary conditions) and an approximately conserved KK-number with Dirichlet boundary conditions. More precisely, in the latter case, c ijk vanishes if i + j + k is even, and falls off as Once we replace the flat space with a slice of AdS, however, additional -10 -JHEP05(2021)096 effects from the bulk will break the KK symmetry, sometimes leaving it approximately conserved as in the flat Dirichlet example, but sometimes not.
The degree and pattern of KK-number violation has an intricate structure. Much of it emerges from the c ijk , through the very interesting properties of the triple Bessel function integrals in eq. (2.13). Specifically, it is convenient to rewrite the integral of interest as a difference of two easier integrals: Except right at threshold, detailed understanding of the I + and I − integrals can be obtained using approximation methods described in section 4. We will use specific cases in the studies of our model presented in section 3. Additional KK-number violation can arise from kinematic constraints. As already noted (see eq. (2.7) and following), certain decays are kinematically forbidden in the single field case for ν > 1/2. These constraints can be more complex in a two field case.
In addition, KK-number might be further violated by interactions at the IR boundary of the space at z = z IR . Specifically, in addition to or as an alternative to the bulk coupling (2.11), we could add an interaction term: This leads to nontrivial interactions if the fields satisfy Neumann boundary conditions at z = z IR . The boundary-localized interaction leads to approximately equal couplings among all the modes, due to (2.10). This contrasts with bulk couplings for small ν i , where the c ijk generally have more structure. A mode's decay width is determined by the available phase space for its potential decays. Recall that the decay width of a scalar where the phase-space function λ PS is defined as We will see cases where near-threshold decays are favored, because approximate KKnumber conservation in the c ijk overcompensates the phase space suppression. This situation leads to near-spherical distributions in decay chains. In other cases this is not so, and decay chains lead to much less spherical events.

Comparison to other high-multiplicity models in the literature
The use of a warped extra dimension to provide a model for a dark or hidden sector is natural following [12,58], and is well-established in the literature [13,[59][60][61][62][63][64]. In this -11 -

JHEP05(2021)096
subsection, we will briefly comment on some related or alternative approaches to modeling high-multiplicity hidden sectors.
Recently, cascade decays in warped dark sectors have been discussed in a series of papers [65][66][67]. While the basic RS framework of these papers is similar to ours, they have focused especially on the regime in which individual KK modes become sufficiently broad that they should be described as a continuum, rather than as narrow resonances. On the other hand, all of our calculations will be done in the regime in which each KK mode is narrow and we can model a cascade decay as a sequence of 1 → 2 decays. This can be done consistently provided that our bulk couplings are sufficiently small, while at the same time not so small that resonances acquire a long lifetime and alter the collider phenomenology. (The long-lived regime may be of independent interest, but its additional complications are beyond the scope of this paper.) In every case that we consider, the couplings can be chosen so that this consistency condition is met. In particular, in every case the width-to-mass ratio of the n th KK mode grows more slowly than it would in a model with decays determined by pure phase space. With pure phase space decays, the n th KK mode has of order n 2 decay modes available to it, with comparable widths. On the other hand, the partial width of a given decay mode scales as in (2.16), i.e., roughly as 1/m n . As a result, the total width of the KK mode scales linearly with the mode number, and the width-to-mass ratio is constant. If we take the bulk couplings to be small but not too small, the width-to-mass ratio of each mode will be small but every unstable KK mode will still decay promptly on collider time scales. The regime of small bulk couplings is where gauge/string duality is well-understood. If the bulk couplings are too small, there is a potential "empty universe problem" in cosmology related to the slow first-order confining phase transition in the dark sector [68][69][70]. However, the details of such a phase transition can be model-dependent (see, e.g., [71]). We assume that potential cosmological problems can be addressed without qualitatively altering the collider phenomenology.
Other recent work [72,73] has studied dark sector particles with simple ansätze for the masses and couplings, e.g., m n = m 0 + n δ ∆m for some positive exponent δ, and couplings depending on factors (m − m i − m j ) r and (1 + |m i − m j |/(∆m)) −s favoring decays to lighter daughter particles or nearby daughter particles, respectively. The model of [72] involves neutral particles χ n decaying toqq χ l , and leads to high-multiplicity event shapes with missing momentum, which are qualitatively similar to events we will discuss. The fully dark decay chains discussed in [73] are even more similar to ours, but are discussed in the context of cosmology rather than collider physics. It is no accident that these decay chains have similar properties, as the ansätze used have been motivated, in part, by models of extra dimensions [74][75][76][77].
Further models for high-multiplicity events have been proposed based on a variety of ideas. SUEP events have been modeled using thermal spectra [15]. Models of black hole production at colliders predict similar spectra [78,79], as do "string balls," lower-mass precursors of black holes [80,81]. Another model with interesting phenomenology, albeit without a well-motivated UV completion, achieves a wide range of couplings within a large ensemble of particles through a random mass matrix [82] (see also [83][84][85]). It would be interesting, in the future, to apply the event isotropy variable to more of these models.

Simulation results
In this section we study event shapes of our toy model for different parameters, by simulating decay cascades of a heavy KK mode with n = n p 1, in its rest frame, to light and stable KK modes. (We will often refer to the KK modes as "hadrons," using the dual viewpoint, but one should keep in mind that these are hidden-sector hadrons, not SM hadrons.) We further decay the hidden-sector stable hadrons (HSH) into two massless particles each, to mimic decays back into the SM. Then we calculate observables from the collection of massless, final-state momenta.
In this simplified model, mass and KK-number are closely related, as we have seen in eq. (2.7). It follows that the degree of violation of KK-number is directly tied to kinematics. Roughly, if KK-number is conserved or lightly violated in a two-body decay, the final state particles tend to be slow in the initial particle's rest frame, while if KK-number is strongly violated, the final state particles are produced with a substantial boost. It is not surprising then that KK-number violation correlates closely with event-shape variables, as we will show in this section.
In extreme limits, it is clear how this should work. Were KK-number precisely conserved in all decays, then every decay would occur at or near threshold, and the final state of the hidden sector cascade would be a collection of slow HSHs. When the HSHs themselves decay to the visible sector, they would produce an array of roughly back-to-back massless particles, produced at random angles. The expected distribution of such particles is roughly spherical. Note, however, that even with dozens of particles, random fluctuations are large and observed events are far from spherical, both to the eye and to event shape variables.
Conversely, in cases with large KK-number violation, the first decay in the cascade produces two relatively light KK states at high boost. Once this occurs, all ensuing decays of the lighter daughters will be highly collimated, and so, independent of the details, two hard jets result. 4 Thus KK-number violating decays early in the cascade leads to a highly non-spherical pattern.
The simulations that we describe in this section interpolate between these extremes. We will demonstrate this using thrust and event isotropy, both of which are sensitive to the features of the events.
These variables are somewhat correlated with a third, namely particle multiplicity. The possible maximum particle multiplicity in our simulations is 2n p . This occurs when KK-number is conserved in every decay, and all hidden hadrons can decay except the n = 1 state, the unique HSH. Then the decay cascade leads to n p HSH's, and to 2n p massless particles once these decay. Violations of KK-number in the cascade, and the existence of multiple HSH's with n > 1, will decrease this number. This tends to increase the event isotropy, since, as noted in eq. (1.5), for a multiplicity k < 192, I sph 192 2/k. (The average isotropy tends to be higher than this estimate, which holds for maximally symmetric events,

JHEP05(2021)096
because of random fluctuations in the angles.) Despite this we will see that isotropy and particle multiplicity are not redundant.
We now present our results, progressing from the most isotropic scenario to the least. For each choice of bulk masses and couplings, we generate 10 4 events starting at KK mode n p = 100 and allow it to cascade into stable hadrons, each of which then decays to a pair of massless particles. In each case, we will see that the degree of KK-number conservation, as reflected in the couplings c ijk and the resulting branching fractions, determines the qualitative properties of the event shapes.
Note that we mainly limit ourselves to small values of ν. This is because scalars with large 5d mass correspond to 4d operators with large scaling dimension (d O = ν + 2), and it is difficult to imagine coupling the SM to them.

Spherical and near-spherical cases
To set a baseline, we begin with the most spherical case in our AdS-based simplified model, a single-field model with ν = 0. This case corresponds to a five-dimensional scalar with mass-squared −2, at the Breitenlohner-Freedman bound, and thus to a dual CFT operator of dimension 2 with a non-zero three-point function. We will compare it a pair of toy models in which KK-number is exactly conserved and spherical events are to be expected. The event shapes for ν = 0 are virtually the same as for the toy models, despite the former's mild KK-number violation and semi-relativistic velocities.
The ν = 0 single field model has a spectrum approximately given by m n ≈ π n − 1 4 . The only HSH is the mode with n = 1; all decays For any single field with a cubic self-interaction and even integer ν, the couplings of its modes satisfy a simple approximate formula. As noted in eq. (2.14), the triple Bessel integral eq. (2.13) can be conveniently written as a difference of integrals. For even ν ≥ 0, I + vanishes due to a factor of 1/Γ(−ν/2) which can be seen in eq. (4.11). Meanwhile I − is approximately given by eq. (4.9). Using eq. (2.9), we find While this is accurate only away from threshold, for ν = 0 the approximation works to within 2% percent except for n 2 + n 3 = n 1 − 1, where the real coupling is smaller in magnitude by up to 5%, and n 2 + n 3 = n 1 , where the difference reaches nearly 30%. From this formula it follows that partial widths behave as λ −3/2 PS , and so decays tend to occur at or very near threshold. This implies that the leading decays for each hadron conserve KK-number. (For instance, the particle with n = n p = 100 has a 77% branching fraction to conserve KK-number, and this varies slowly with n p .) Even those decays that violate KK-number do so by small amounts, and in the end the average HSH multiplicity at the end of the cascade is reduced only to 93 from its maximum of 100; this is shown later in figure 5. The decays of the HSHs produce nearly 200 massless particles with an energy distribution, shown in figure 3(c); note that m 1 ≈ 2.40 ≈ m 100 /130, and the distribution peaks at about m 1 /2, with a tail up to ∼ 2m 1 . Thus velocities of the HSHs tend to be only semi-relativistic, and the angular distribution of their massless daughters is largely random. Our first toy model has the same masses and KK-number-conserving couplings as the ν = 0 case just described, but we set all KK-number violating couplings to zero by hand: We will call this the KK-conserving (KKC) ν = 0 model. Its final state consists of exactly 200 massless particles, with an energy distribution slightly narrower than the full ν = 0 model, as the latter has KK-number-violating decays with more kinetic energy. We expect it to have slightly more spherical events. The second toy model (the "flat case") is a single field on Minkowski space times an interval, M 4 ×S 1 /Z 2 , with Neumann boundary conditions on the field. As noted earlier, the spectrum has m n ∝ n exactly, and KK-number is conserved. Strictly speaking all particles are marginally stable, but we imagine deforming the model infinitesimally so that all decays can occur. The final state from an initial heavy hadron with quantum number n p consists of n p HSHs with n = 1, all at rest. The decay to SM massless particles then produces events with exactly 200 massless particles in back-to-back pairs, distributed randomly in angle. Each particle has energy exactly m 1 /2 = π 2 .

JHEP05(2021)096
Now we compare the event shapes for these three cases. In figure 3 we show the distributions in energy, event isotropy, and thrust for the particles in the final state. For consistency in range, we plot the scaled thrustT such that all variables have a range of [0, 1] with 0 being the most isotropic and 1 being the least. All three of these examples are very similar as seen by event-shape variables. The most notable differences are percent-level shifts in event shapes, in tails that arise from small numbers of somewhat less isotropic events. We may therefore treat any one of them as a benchmark against which to compare other cases.
All three cases have event isotropy that peaks in the range 0.15-0.20. From eq. (1.5), maximally isotropic events with 192 particles U sph 192 would have Naively we might have expected the flat case, with 200 particles of equal energy and random angles, to approximate this value. However, the random fluctuations in angle (but not in energy, which remains m 1 /2 for each particle) lead to a significantly higher event isotropy, closer to 0.16; we do not know a method to compute this number without simulation. 5 Despite the wider energy distribution of the ν = 0 and KKC cases, their event isotropies are quite similar to the flat case. We will explore the causes of this in [54].

Single field, general ν
Next, still studying a single scalar field with Dirichlet boundary conditions, we consider other values of ν. Each has a different degree of KK-number violation. Although the amount of KK-number violation is still relatively small, and the decays are still mostly close to threshold, the effects are large enough to observably shift event shapes relative to the ν = 0 benchmark. As noted in eq. (2.14), the integral eq. (2.13) can be written in terms of a difference of two integrals I + and I − . Substituting ν i = ν into eq. (4.13) and using the Euler reflection formula, one finds where H is a function whose dependence on ν and the m i is subleading compared to the terms shown explicitly. 6 When ν is an even integer, I + vanishes, so the couplings are determined entirely by I − . Since, from eq. (4.9), I − ∼ λ −1 PS near but below threshold, the branching fractions decrease as ∼ λ −3/2 PS , strongly suppressing KK-number violation. For ν = 0, where there 5 For KK-number conserving scenarios, events become more spherical at high np, as shown in appendix B, and the theoretical limit is reached at large np. 6 The full integral I in eq. (2.14) is elementary at ν = 1 2 , and the λ −1/2 PS factor is easily seen there.

JHEP05(2021)096
are no other sources of KK-number violation, this gives quasi-spherical events. For general ν, however, I + ∼ λ −1/2 PS falls off more slowly away from threshold than I − . Although for ν very large this does not matter, because the power of the mass ratio in eq. (3.5) decreases rapidly with ν and tends to disfavor decays to light hadrons, for ν ∼ 1 one finds |I + | ∼ |I − |, and so KK-number violation in the couplings is considerably larger than for ν ≈ 0, 2, 4, . . . .
The dominant decays remain KK-number conserving for ν 0.1. For the ν = 0.15 case below, KK-number violation in the couplings is already significant, and reduces the multiplicity of HSHs well below the ν = 0 benchmark, as we will see shortly.
For ν ≥ 0.5, a new effect reduces multiplicity further: m n < m n + m n−n , so KKnumber conserving decays all become kinematically forbidden, and the minimal amount of KK-number-violation per decay is > 0, even for even integer ν. This reduces the number of typical decays in the cascade and the number of hadrons at the end of the cascade. However, even though KK-number conservation is forbidden, decays with large KK-number violation are still somewhat suppressed, and so the leading decays are minimally KKnumber-violating -that is, they have the smallest amount of violation consistent with kinematic constraints. These decays are generally the ones closest to threshold.
When decays can typically only violate KK-number by a small amount (per decay), they remain near kinematic threshold, so boosted hadrons and ensuing jetty structures in the events do not arise. Nevertheless, I sph 192 increases. In part, this is due to a decrease in particle multiplicity. The total amount of KK-number ∆ KK,tot lost in the cascade (equal to the sum over decays of the KK violation ∆ KK,s in each decay s) is given by the parent mode number n p at the beginning of the cascade minus the sum over KK-numbers n i of the N HSH HSHs in the final state, where S is the total number of decays in the cascade. The reduced final-state hadron multiplicity N HSH increases the minimum value of I sph 192 by a factor of order 100/N HSH , from eq. (1.5), even if the HSHs are rarely boosted and their final decay products are quasi-isotropically distributed.
If KK-number is exactly conserved, S = n p − 1 and N HSH = n p . More generally N HSH ≤ n p − ∆ KK,tot , with the equality holding only if the only HSH has n = 1. This is true for ν < 1 2 , and in particular for the case ν = 0.15 that we show below. For larger ν, there are HSH's with n > 1, so N HSH is even smaller; and on top of this, decays with ∆ KK,s = 0 are forbidden, so ∆ KK,tot is of the same order as S and n p , leading to a substantial reduction in N HSH . For the ν = 0.75 case we show below, S ∼ n p /2 and N HSH < n p /2, leading to a reduction of the multiplicity by more than half compared to the ν = 0 case, and a corresponding substantial increase in I sph 192 . The branching ratios for the n p = 100 mode into daughter modes n 1 , n 2 with ν = {0, 0.15, 0.75} are shown in figure 4. For convenience we will refer to these plots throughout the paper as "branching fraction triangles." For ν = 0, the only HSH is the  allowed, but the probability of conservation in the decay of heavy modes is only ∼ 30%. Decays with KK-number violation ∆ KK = 1 occur with comparable probability (∼ 30%). Technically, this is due to a cancellation in the integral for c ijk ; the integrals I + and I − are similar in magnitude for ν ∼ 0.1 − 1.9 and interfere destructively (constructively) when i + j + k is even (odd). The checkerboard pattern in figure 4 (b,c) arises from this effect. The pattern of branching ratios for ν = 0.75 is similar to that of ν = 0.15. However, ∆ KK = 0 decays are kinematically forbidden, as can be seen by careful examination of the upper right edge of the triangle.
While figure 4 applies for n p = 100, it illustrates qualitative features that apply, at fixed ν, for smaller values of n, and thus for the whole cascade. This is because the integrals I + and I − have relatively simple behavior under changes of n; see appendix B for some discussion. We present results for the event shapes of the cascades in figure 5, where we show distributions of the multiplicity and energy of the massless particles in the sample, along with the sample's event isotropy and thrust distributions. KK-number violation reduces the multiplicity for ν = 0.15 relative to the ν = 0 benchmark. The ν = 0.75 case has even smaller multiplicity and a wider energy distribution, due to additional sources of KKnumber violation and the fact that the n = 2 state is also an HSH. Meanwhile the event isotropy increases (the events become less spherical) for ν = 0.15 and even more so for ν = 0.75. Interestingly, this pattern does not apply for thrust; the thrust distributions for ν = 0.15 and 0.75 are very similar. This suggests that event isotropy and thrust are sensitive to different event shape characteristics and are not redundant variables. Further study of this issue will be presented in [54].
Comparing figure 5 (a) and (c), one might wonder if event isotropy and particle multiplicity are redundant variables, especially considering the natural correlation between them that we described earlier. Although the correlation seems particularly strong in these examples, other simulations shown below clearly demonstrate these two variables are independent.
We illustrate some characteristic events for each of these samples in figure 6. An event of average isotropy in the most isotropic sample (ν = 0) does not show any clear boost We end our discussion of the single field case with a few comments on large values of ν. With increasing ν but fixed n p , I + is suppressed by a factor of m 2 m 3 /m 2 1 ν , so I − dominates in most decay channels. When this happens, the couplings conserve KKnumber similarly to eq. (3.1) even when ν is not an even integer. However, this fact is irrelevant for n p = 100 because the kinematic constraints imposed by the mass spectrum require large violation of KK-number in each decay, which grows with ν. We noted earlier that KK-number conservation becomes kinematically disallowed for ν > 1/2. More generally, decays with KK-number violation ∆ KK begin to be kinematically constrained when ν 2∆ KK + 1/2, starting with the most symmetric decays, and once ν reaches the values 0.5, 2.92, 5.56, 8.30, 11.1, . . . , all decays with ∆ KK = 0, 1, 2, 3, 4, . . . are forbidden. Moreover, there are roughly 1 + 1 2 (ν + 3 2 ) HSHs in the spectrum. All of these effects drastically reduce multiplicity and widen energy distributions, further increasing the event isotropy and thrust, so we do not expect highly spherical distributions to be common in this regime.
Because there are two towers of hadrons with different mass spectra, the correlation between KK-number and mass is not as simple as in the single-field case. However, these complications are relatively unimportant compared to the dramatic change in the pattern of couplings c ijk . In the single field case, we have seen that couplings near threshold, with zero or minimal KK-number violation, are strongly enhanced, because I − ∼ λ −1 PS (away from threshold) and I + ∼ λ −1/2 PS . This leads to events that, to a greater or lesser extent depending on ν, tend to be quasi-spherical; jetty events are rare. But this behavior does not extend to general ν i .

JHEP05(2021)096
For the two field case, once ν for a decaying particle is significantly larger than the sum of the ν i for its daughters, KK-number violation becomes large. We can see this by examining I + (which typically is much greater than I − in this regime) using eq. (4.11) and the paragraph following it; see also section 4.3.2. If, without loss of generality, we take the decaying particle to be from field Φ 1 and its daughters from fields Φ 2 , Φ 3 , and set ∆ν ≡ ν 1 − ν 2 − ν 3 , then when ∆ν = 2k, where k is any positive integer, I + goes to a constant at threshold. Consequently branching fractions are suppressed near threshold by λ +1/2 PS , and instead peak elsewhere in the branching fraction triangle. Moreover, I + has k − 1 lines of zeroes, and so the branching fraction triangle has k plateaus separated by valleys. Some of these plateaus have large KK-number violation. When ∆ν = 2k these plateaus survive, but are supplemented by a return of the λ −1/2 PS behavior near threshold. Although this near-threshold enhancement favors KK-number-conserving decays, the large KK-number violating decays in the more distant plateaus often remain dominant, making jetty events common. 7 Focusing now on ν 2 = ν 3 = 0, we will illustrate the behavior just described in the cases ν 1 = 2, 3, 4, whose branching fractions for n p = 100 are shown in figure 7. For ν 1 = 2, using eq. (2.9), eq. (2.13), eq. (4.11), and the remark below Eq. (4.11) that F 4 → 1 in this case, we find (to a very good approximation) The branching fractions are then proportional to m 2 m 3 λ 1/2 PS , so they vanish at all boundaries of the branching fraction triangle and are broadly distributed, as seen in figure 7(a), peaking ∼ 10% below threshold. For ν 1 = 4, eq. (4.12) gives which creates a zero between the threshold region at upper right and the m 2 , m 3 → 0 corner at lower left. The branching fraction triangle then has two plateaus, one far from threshold and one nearby, as seen in figure 7(c). The probability for a particle to decay via either plateau is comparable. Finally, for ν 1 = 3, the zero seen for ν 1 = 4 is still present, closer to threshold, but in addition there is enhancement right near threshold. Despite this enhancement, there are so many decay paths in the plateau that the total probability to decay at or near threshold is only ∼ 1/4, and so large KK-number violation is the norm. A cascade in the two field case involves both φ 1,i → φ 2,j φ 2,k and φ 2,i → φ 1,j φ 2,k decays. To compute the branching fractions for the latter, we need to exchange ν 1 and ν 2 (but not 7 In this discussion we have neglected I−. It can be seen from the formulas of section 4 that |I+| |I−| for most decays, but there is a subtlety near threshold, where our approximation I− ∼ λ −1 PS seems to blow up faster than I+ does. This effect is merely due to our approximation eq. (4.9), however, which is not valid at threshold. Instead, when m1 − m2 − m3 → 0 for some decay, which requires tuning of the νi, both I+ and I− ∼ λ −1/2 PS extremely close to threshold, as noted in eq. (4.10), and in fact the original integral in eq. (2.13) is always finite there. (See also appendix B.) In practice, then, our discussion here of I+ captures all the important features of the branching fraction triangles, except for the precise details in the near-threshold region.

JHEP05(2021)096
ν 3 ) in our analytic formulas. For reasons similar to those leading to eq. (3.5), I + vanishes for ν 1 = 2, 4, giving the nearly KK-number-conserving result eq. (3.1), but is important for ν = 3, leading to slightly higher KK-number violation. Thus it is no accident that the branching fraction triangles for these decays, shown in figure 8, resemble those of the single field case, figure 4. However, the KK-number violation in φ 2,i decays is subleading compared to the much larger KK-number violation that can occur in φ 1,i decays, and its details do not much impact the results.
As was true also for figure 4, the qualitative features of the branching fraction triangles are present also for smaller values of n p , and thus apply for the whole cascade. See appendix B for further discussion of the n p dependence.
Particle multiplicities are affected not only by these KK-number-violating processes but also by the increasing number of HSHs. For ν 2 = ν 3 = 0, and ν 1 1.75, the decay φ 1,1 → φ 2,1 + φ 2,1 is always open, so only states of Φ 2 are stable against decays within the hidden sector. The mode φ 2,1 , as the lightest mode in the hidden sector, is of course an HSH, while φ 2,2 is stable for ν 1 ≥ 0.5, φ 2,3 is stable for ν 1 2.90, and φ 2,4 is stable for ν 1 5.53. But the effect on event shapes of reduced multiplicity is subleading compared to jet creation through boosts of daughter particles, as we will now see.
Our results are shown in figure 9, based again on cascades starting at the 100 th KK mode of the field Φ 1 , with 10 4 events per sample. Just as the branching fractions are much more widely distributed in figure 7 than in figure 4, all the distributions in event isotropy and thrust are much wider than for the single field case in figure 5. As is evident by comparing these figures with figure 7, the plateaus and valleys in the branching fraction triangles of the ν 1 field lead directly to structure in the event shape variables.
It is easy to see why this is so. The early stages in the cascade are most important, because decays far from threshold early in the cascade create boosted particles with a substantial fraction of the event energy, and their collimated decay products lead to hard jets, while reducing the overall multiplicity of hadrons. Such events will differ strongly from the near-spherical events we saw in the single-field examples. If instead the initial decays nearly conserve KK-number and create several slow heavy hadrons, KK-number violation in their decays will still lead to jets, but these will have a much smaller fraction of the event's energy. Thus, depending on the details of a particular cascade decay, an event may present a small number of hard jets at one extreme, or a small number of soft jets on top of a quasi-spherical background at the other. As the probability for far-from-threshold decays increases, so will the fraction of events with larger thrust and event isotropy. This is seen most dramatically for ν 1 = 4, where the bimodal feature in the branching fractions, figure 7(c), with one plateau that weakly violates KK-number and another that violates KK-number significantly, leads to bimodal event shapes at small and large event isotropy and thrust respectively. These correspond to two classes of events, one relatively spherical and the other relatively jetty. Note that the lobe at larger event isotropy is comparable to the event isotropy distribution for threshold tt events, shown in figure 1, and actually peaks at higher I sph 192 , though it is not as jetty as a qq sample. Comparing ν 1 = 2 and ν 1 = 3, one sees the latter's ensemble of events is less spherical on average. Despite the near-threshold enhancement for ν 1 = 3, the majority of decays occur in the plateau far from threshold. This plateau is located further from threshold than that of ν 1 = 2, so the most common decays for ν 1 = 3 have larger boost on average than those for ν 1 = 2, and this leads to harder jets, moving the isotropy and thrust distributions to larger values. The near-threshold decays partly compensate for this effect, making the event-shape variable distributions particularly broad. The multiplicities of massless particles are of order 80 for ν 1 = 2 and (partly because of the additional HSH) only 50 for ν 1 = 3, 4. Although this was not evident in the single-field cases, here one can see clearly that the event shape variables are not perfectly correlated with multiplicity. For instance the multiplicity distribution for ν 1 = 4 does not show the obvious two-lobe structure seen in the event shape variables. We will return to this issue briefly in section 3.5. It is also noteworthy that the isotropy distribution is narrower than the thrust distribution, a hint of imperfect correlation which we will study further in [54]. Turning our attention to the energy distributions of figure 9, it is interesting to note that all the spectra (including the reference case ν = 0) are peaked at the same value. The peak is at half the lightest mass in the cascade: m ν=0 1 /2 = 1.2. This tells us that for all cascades there are many soft particles produced in the decay of nonrelativistic n = 1, ν = 0 particles. The width of the distribution is much wider for ν 1 = 2, ν 2 = 0 than for ν = 0, and even wider for ν 2 = 3, 4. This tail is due to the fact that modes heavier than the lightest mode are stable in the two field cases, whereas for ν = 0, only the lightest mode is stable. Note that these observations also apply to all the single field cases shown in figure 5.
Visualizations of characteristic events from each sample are shown in figure 10. The effects of the unsuppressed KK violating decays are evident.
Reference Case: Reference case: Reference case: Figure 9. Same as figure 5 but for the two field scenario with ν 2 = ν 3 = 0, and ν 1 = 2 (blue), ν 1 = 3 (green), and ν 1 = 4 (red). Note that average multiplicity is lower and event isotropy is higher than the single field cases in figure 5. The appearance of plateaus in the branching ratios can greatly increase the amount of KK violation per event.

Cascades with boundary couplings
In the scenarios discussed above, we have assumed that fields interact in the bulk of the extra dimension, as in (2.11), and we have taken the boundary conditions to be Dirichlet. A very different phenomenology arises if we assume that the fields have Neumann boundary conditions, and add the interaction term (2.15) localized on the IR boundary of the space. As indicated in (2.10), the boundary values of the different KK modes are approximately equal, and so decays are governed approximately by the phase space of the final state. The branching ratios in this case are illustrated in figure 11, which also shows the (similar) branching ratios that would result from pure phase space factors. Of course, one can also consider a model containing both the bulk coupling c (2.11) and the boundary couplingc (2.15), obtaining physics that interpolates between the two. Below we will show some results from such a case withc = 0.015c, which turns out to produce event-shape distributions roughly midway between those of the pure bulk and pure boundary cases.
In figure 12, we illustrate how boundary couplings affect the event shape. The blue curve shows the case of pure bulk couplings with Neumann boundary conditions, for a single self-interacting field with ν = 0.3; as in figure 5, this leads to approximately spherical events with small values of I sph 192 and thrust. 8 The orange curve shows events with boundary decays, which are substantially less isotropic. They also have a broader distribution of both event isotropy and thrust, although the event isotropy distribution is more peaked than the thrust distribution. Comparing to figure 1, we see that the typical event with boundary decays is more isotropic than a QCD dijet event, but has similar isotropy to a near-threshold tt event. Finally, we show the case with a mixed bulk/boundary coupling,c = 0.015c, in green. As expected, this interpolates between the bulk and boundary cases. It does so by broadening the distribution, rather than producing a narrow peak in between the two cases. Visualizations of typical events are shown in figure 13. figure 12 also shows the typical multiplicity and energy of the individual final-state massless particles in the events. We see that the boundary cascades have much lower multiplicity, because the decays more often go directly to lighter daughters, so it takes fewer steps to reach the HSHs at the bottom of the cascade. Consequently, the individual particles also have more energy than they would for bulk couplings. This figure suggests that the event isotropy may be highly correlated with the particle multiplicity. Although this is true when we start all cascades with the same initial KK-number, it is not the case in general: boundary cascades remain much less isotropic even if they begin with a much larger choice of n p . Thus event isotropy captures different information from particle multiplicity, or even from pairs of observables like particle multiplicity and thrust. We will discuss this more in the companion paper [54].

Evolution of multiplicity, event isotropy in cascade
To investigate further the correlations between multiplicity and isotropy, it is interesting to see how these evolve through the decay cascade. We do the following exercise.  Figure 11. The branching ratios for ν = 0.3 of the 100th KK mode into all kinematically allowed two-body final states, in scenarios with Neumann boundary conditions. The cases are (a) bulk couplings; (b) boundary-localized coupling; and (c) pure phase space. In (d), we show a 1d plot of the branching ratios along the n 1 = n 2 line. For the boundary-localized case (b), in contrast to figure 4 and the bulk case, the decays significantly populate the full triangle, and favor daughter particles with substantial momentum. The distribution is similar to that of pure phase space (c). step in the cascade, we take all the hadrons present at that step (independent of whether they are stable against hadronic decays in following steps) and artificially force them to decay to massless particles. At the initial step of the cascade, with just the n p mode at rest, every event has a pencil-dijet with I sph 192 =T = 1, while at the end of the cascade we obtain the samples studied above. In between, the multiplicity in each event gradually increases and I sph 192 decreases. In figure 14, we plot the average multiplicity and the average isotropy (averaged over 10 4 events) at each step in the cascade. We also show, as a dashed line, the theoretical lower limit on I sph 192 for the corresponding multiplicity; see eq. (1.5). If multiplicity and event For ease of comparison with other figures, we also show the ν = 0 Dirichlet case (originally displayed in figure 3). We see that the cascade with a boundary coupling produces less isotropic events (larger I sph 192 ). Turning on both couplings interpolates between the bulk and boundary case, with a broad distribution of I sph 192 .
isotropy were perfectly correlated, then all of the cascades would lie on the same curve, even though at any given step, and at the end of the cascade, they would sit at different values. Instead, we see that the cascades for different choices of ν i can give different curves. This is additional evidence that isotropy measures more than multiplicity. We will explore the independence and correlation of multiplicity, isotropy, thrust and other event shapes in [54].

Analytical estimates for couplings
In section 3, we saw that different parameter choices lead to qualitatively different patterns of couplings c ijk and branching fractions, and from there to qualitatively different event shapes. In this section, we provide analytic calculations of the overlap integrals that determine the couplings, in order to substantiate and explain the results presented above. In particular, an inverse dependence on the phase-space function λ PS will be manifest in our results, explaining the cases in which we have observed a preference for near-threshold decays. We will also understand the oscillatory behavior that gives rise to the observed plateaus in the two-field case.

Basic ingredients and general strategy
Our goal is to gain an analytic understanding of integrals of the form  figure 4) along which we have done the calculation. At left, we take near threshold decays (right-hand edge of the triangle); at right, we take decays along a slice through the middle of the triangle. For near-threshold decays, we see that I + and I − are comparable. Away from threshold, I + dominates, with I − contributing a small oscillatory pattern.
We choose a convention where particle #1 is the heaviest particle, i.e., we assume without loss of generality that m 1 ≥ m 2 + m 3 . We can separate our integral into two pieces, (4. 2) The reason for doing this is that when all of the masses (in units of the IR brane scale) are sufficiently large, i.e., when m i ν 2 i , we can approximate the Bessel functions in the integrand of I − (ν i , m i ) by their large-argument asymptotic expansions. This makes approximating I − (ν i , m i ) into an analytically tractable problem. On the other hand, the integral I + (ν i , m i ) over the whole positive real axis is known analytically. By combining the exact analytic answer for I + and the approximate answer for I − , we obtain an analytic approximation to the I(ν i , m i ) and hence to the couplings among Kaluza-Klein modes.
In figure 15, we show examples of how I breaks down into contributions from I + and −I − , along two slices of the branching ratio triangles in the single-field case with ν = 0.75. Typically, I + dominates, and the integral is suppressed for decays to light modes. Near threshold, I − gives a contribution comparable to that of I + . As we will see below, in certain special cases (like ν = 0), the integral I + is identically zero, but the behavior seen in the plot is representative of more general ν values.
Our goal in the remainder of this section is to provide some analytic insight into the behavior of the integrals I + and I − . We will first discuss an analytic approximation to I − (ν i , m i ), and compare it to numerical results. Then we will present an exact analytic formula for I + (ν i , m i ), and comment on a special case to elucidate the "plateau" structure observed in figure 7.

Approximating the integral I − (ν i , m i )
We apply the large-argument asymptotic approximation of the Bessel function, (2.6), to estimate I − (ν i , m i ): We can rewrite the product of three cosines as a sum of four cosines by repeated use of the identity 2 cos a cos b = cos(a + b) + cos(a − b). Given signs σ, σ ∈ {+1, −1}, we define In labels, we will suppress the "1" and write the σ's as + or −; for example, we denote m 1 + m 2 − m 3 by m +− . Then (4.3) is equivalent to (4.5) This integral can be performed analytically in terms of the Fresnel cosine and sine integrals. Our conventions for these are specified in appendix A.1. We obtain: (4.6) This is already a useful approximation to I − (ν i , m i ). We can go further by noting that, provided the masses are all large and that we do not consider decays too close to threshold, we can exploit the large-argument asymptotics of the Fresnel integrals, (A.2). In this case, we find that the answer depends on sin and cos of m σσ , so we make use of the mass eigenvalue estimates in (2.7). We will provide the estimate in the case of Dirichlet boundary conditions.

I − (ν i , m i ) estimate for Dirichlet boundary conditions
Keeping the first subleading, O(1/z), term in the Fresnel integral asymptotics, the expression in brackets in (4.6) becomes Using the sine sum-of-angles identity, then using (2.7) to replace the masses in the argument of the sine by an approximate expression in terms of the KK mode numbers n i , we obtain -31 -

JHEP05(2021)096
a simple approximate formula for the Dirichlet case: where n i is the integer mode number of the corresponding KK mode. In particular, there is a term that scales as 1/(m 1 − m 2 − m 3 ) that dominates when the phase space is relatively small (but not so small that our approximations break down). We can write this term using the phase space factor (2.17) , (4.9) which we used in (3.1). The large-argument expansion of the Fresnel functions breaks down at m σσ 1, very close to threshold. This case generally only arises if we tune the values of ν so that the offsets in the Bessel function zeros align to allow for m 3 ≈ m 1 + m 2 , so we do not expect that it is generally relevant, but we discuss it for completeness. In this case, the sum is dominated by the small-argument expansion of the Fresnel functions for the σ, σ = −1 case, . (4.10) Thus, very close to threshold, we expect that the divergence is ameliorated to λ −1/2 PS , except in cases where ν 1 − ν 2 − ν 3 is an even integer, when this term has coefficient zero and subleading terms dominate. As we will see below, this is also the near-threshold behavior of the integral I + , so we predict that I − does not parametrically dominate over I + in the small phase-space region. In fact the original integral is finite at threshold, so the singular behavior of I + and I − must cancel there.
Unlike the Dirichlet case, note that the large-argument approximation for the Fresnel function would have given zero for the Neumann case, because the different constant term in (2.7) removes the π 2 (1 + σ + σ ) term in the argument of the sine and leads to zero. To obtain a similar approximation in the Neumann case, we must keep subleading terms in the various approximations we have made. We will not do so here.
We compare a numerical computation of the integral I − with the two approximations I − , works very well away from threshold but deviates close to the threshold, leading to the visible dotted curves. As in figure 15, the inset triangles illustrate the slice through the branching ratio triangle that is plotted. We show the case ν = 0.3 because it shows a larger discrepancy, and thus a more visible dotted curve, than the case ν = 0.75.

General formula
In the case of I − , we used the asymptotic expansion of the Bessel function to obtain an analytic approximation. We can do better with I + : the integral is analytically known (eq. (7.1) of [86], in the special case λ = 2) to be: (4.11) The function F 4 is known as an Appell function; it is a two-variable generalization of a hypergeometric function. This integral has previously appeared in the physics literature on 3-point correlators in momentum space in conformal field theories [87][88][89] and de Sitter space [90,91]. For convenience, we include its definition in appendix A.2. The behavior of the Appell function leads to one of the important qualitative features we have observed in our numerical results: the existence of plateaus of large branching fractions separated by valleys of suppressed branching fractions. It is manifest from the series definition of the Appell function that it becomes a polynomial when its first argument is a non-positive integer. When ν 1 = ν 2 + ν 3 + 2k, for k a positive integer, the polynomial has degree k −1, and has k −1 curves of zeroes. Numerical results, supported by incomplete analytic arguments, indicate that these zeros always lie in the physical region m 2 + m 3 < m 1 . 9 Two special cases of interest to our two field case are k = 1, for which the Appell 9 Along the line m3 = 0, this follows from the fact that F4 is simply 2F1, discussed in appendix A. (4.12) The zeros of the Appell functions form the valleys between plateaus observed in figure 7. Even in cases where the Appell function is not a polynomial, we expect that ν 1 −ν 2 −ν 3 2 , when positive, approximately counts the number of plateaus.
The expression eq. (4.11) can always be reduced to an expression in terms of ordinary hypergeometric functions [93][94][95][96] . (4.13) The variables X and Y are defined so that ( We comment on some details of the reduction from F 4 to 2 F 1 in appendix A.3. We can extract from (4.13) a few simple, general points: • The integral I + has a phase-space factor λ 1/2 PS (m 2 1 , m 2 2 , m 2 3 ) in the denominator. Thus, precisely when the general decay width formula (2.16) would lead one to naively expect a decay to be rare, the coupling squared enhances it via an inverse dependence of the decay width on the phase space. This accounts for the general tendency of the 5d models to feature many near-threshold decays.
• The hypergeometric functions have the property lim z→0 2 F 1 (a, b, c; z) = 1, so they are unsuppressed when X or Y is small. However, the prefactors (m 2 /m 1 ) ν 2 (m 3 /m 1 ) ν 3 indicate that, in general, the decay rates to light daughters are suppressed. Even in the special case ν 2 = ν 3 = 0 when this suppression is absent, the normalization factor N (ν) n is smaller for light daughters than heavy ones, as indicated in (2.9). Thus, in general, we expect that decays to light daughters are rare, as confirmed by the numerical results in figure 7, for example.
• Recall that Γ(x) has poles whenever x is a nonpositive integer. As a result, whenever ν 2 + ν 3 − ν 1 is a positive even integer, the integral I + (ν i , m i ) vanishes. In the special together with Theorem 3.2.i of ref. [92]. In the special case ν2 = ν3, it can also be proven along the line m2 = m3, using the same theorem together with eq. (4.13) below. However, we lack a completely general proof.

JHEP05(2021)096
case ν 1 = 0, there are also zeros in the denominator, but taking the limit from nonzero ν 1 's shows that I + vanishes in this case as well (as we saw in figure 15, where all ν's were zero).
• We have X, Y ∈ (0, 1), so that the hypergeometric functions are nonsingular in the physical region except perhaps near threshold. In the near threshold region m 2 + m 3 → m 1 , we have X → m 2 /m 1 and Y → m 3 /m 1 . The hypergeometric functions may be singular when X → 1 or Y → 1. For concreteness, consider the case X → 1. This requires that m 3 → 0, m 2 → m 1 . In this case, the 2 F 1 functions of X diverge as (1 − X) −ν 3 . However, this is compensated by the prefactor (m 3 /m 1 ) ν 3 . As a result, the prefactor of λ −1/2 PS is the only source of singularities at the boundary of the physical region.

4.3.2
Plateau structure in the special case ν 2 = ν 3 = 1/2, ν 1 ∈ Z We already noted that F 4 becomes a polynomial when ν 1 −ν 2 −ν 3 is a positive even integer; in this case, the hypergeometric functions in (4.13) are also simply polynomials in X and Y . We will now present a special case for which we can give a straightforward derivation of the integral I + in terms of Chebyshev polynomials. This is one of the simplest cases in which the existence of plateaus and valleys can be deduced.
Spherical Bessel functions have simple expressions in terms of trigonometric functions. In particular, for ν = 1/2, we have the identity (4.14) In the special case ν 2 = ν 3 = 1/2, this reduces our overlap integral to (4.15) The infinite integral ∞ 0 dt J ν (at) cos(bt) is known [97, section 13.42]. Here, we will provide a clear derivation in the special case when ν is an integer. The result is 0 if ν is odd, whereas if ν is even and b < a, it is given by [98,99]. The link between Chebyshev polynomials, defined by T n (cos φ) = cos(nφ), and Bessel functions of integer index arises from the decomposition of plane waves in cylindrical coordinates (ρ, z, φ): (4.16) One could take this to be a definition of the integer-index Bessel functions via a generating function. Sending i → −i, we learn that J −m (ρ) = (−1) m J m (ρ). By integrating both sides of (4.16) against cos(nφ) from −π to π and exploiting orthogonality, we obtain the integral representation J n (ρ) = i −n π π 0 dφ e iρ cos φ cos(nφ), where we have used that the cosine is even to rewrite the integral from 0 to π instead of −π to π. In this expression we see that taking a derivative with respect to ρ brings down a factor of i cos φ. Using the fact that cos(nφ) = T n (cos φ), we then obtain the identity This will allow us to easily compute the integral ∞ 0 dt J n (at) cos(bt) for arbitrary integer n once we know the integral for the special case n = 0.
For n = 0, we directly use the representation (4.17) of J 0 : One way to see this is to view the integral over φ as a contour integral, which picks up poles where cos φ = ±b/a which are in the domain of integration (−1 ≤ cos φ ≤ 1) when b < a but not otherwise. Next, we can obtain the result for general J n by using (4.18) and then integrating by parts to move the derivatives onto the cos(bt) factor. If n is even, T n (x) contains only even terms and the derivatives produce a cos(bt) factor; if n is odd, a similar argument leads to a sin(bt) factor. Every two derivatives acting on a cos(bt) factor will multiply by −b 2 = (ib) 2 , effectively absorbing an extra factor of ib into the argument of the polynomial. Hence: agreeing with the results in the literature [98,99]. The polynomial behavior of T 2k (b/a) leads to several zeros as a function of b/a, and hence to "plateau" structure in the Bessel overlap integrals like that we have previously observed (for different choices of ν) in figure 7. More generally, the plateau structure arises due to similar oscillatory behavior in the hypergeometric functions in (4.13).

Conclusions
One of the challenges facing LHC studies is that new physics with unusual signatures might be able to escape trigger strategies or hide in the large data sets. It is important therefore to consider these unusual signatures carefully, especially those that rarely appear in classic BSM models and are difficult or impossible to calculate with confidence. Among these signatures are complex high-multiplicity final states. Events with a small number of QCD-like jets are well-studied, and various approximately spherical high-multiplicity signals have also been considered, but little is known about signatures that, in some sense, lie between these extremes.

JHEP05(2021)096
It is well-known that large 't Hooft coupling gauge theories can produce spherical events [13,32,33], and that RS-like 5d models with cascade decays of Kaluza-Klein modes can produce approximately spherical events [52]. Here we have shown that, in fact, 5d simplified models with tunable parameters (including a small number of bulk fields with bulk or boundary interactions) can produce a wide range of event shapes in cascade decays of their heavy states. These 5d simplified models are well-suited to serve as templates when designing collider searches for unusual events that might be hiding in samples of events with high jet multiplicity.
Specifically, we saw that a key determinant in the cascade decays was the degree to which KK-number is violated, which in turn determines how close to threshold are the majority of decays. With some choices of bulk parameters, KK-number is approximately conserved, leading to quasi-spherical event samples; see figure 5. For other choices of bulk parameters, or with the addition of boundary interactions, KK-number can be strongly violated, making events with a few hard jets commonplace. In the latter cases, one can find samples as jetty as threshold tt events, with many individual events as jetty as SM qq events; compare figure 9 or figure 12 to figure 1. In demonstrating this, we have relied not only on thrust, a classic event shape variable, but also on event isotropy, a newly-introduced variable which appears well-suited to this purpose.
One of the main results of this paper is an approximate analytic understanding of the Bessel function overlap integrals eq. (4.1) that determine the couplings among different KK modes in the 5d simplified model, which are the most important source of the KK-number violation. We are not aware of any previous detailed studies of this definite integral. In particular, we have seen that, aside from decays very near threshold, the integral over the finite fifth dimension is generically well approximated by the integral I + over an infinite interval; see eq. (4.2). This integral has an analytic expression, (4.13), in terms of hypergeometric functions. This integral has a factor of the phase space function λ 1/2 PS , defined in (2.17), in the denominator. As a result, the naive expectation that decays far from threshold are favored due to the larger available phase space is precisely inverted within the context of these extra-dimensional models. Decays with small phase space are often favored. Similar results hold for the complementary integral I − , which dominates in special cases, and for which we have provided an approximate analytic understanding. However, there are also cases where I + exhibits quasi-polynomial behavior with multiple zeroes, and is consequently enhanced far from the threshold region and/or suppressed at the threshold region. Then decays near threshold are no longer dominant and the cascades are more likely to produce kinematic jets.
In a companion paper [54], we will explore the event shapes that arise from our 5d simplified models in more depth. In particular, we will illustrate, using both simulation and analytic estimates, that event isotropy provides an important complementary probe, capturing aspects of event shapes that are distinct from those captured by thrust, the eigenvalues of the sphericity tensor, or jet multiplicities.
A remaining task is to connect our 5d simplified models with the Standard Model. In this paper, we have started our events with a single heavy KK mode which then cascades into many daughter particles. In order to use event shape observables based on massless JHEP05(2021)096 momenta, we have assumed that all of the final-state daughter particles decay into two massless particles. However, we have stopped short of a full model for the interaction with the SM, which would allow for the production of many different modes and for a more general set of decays. More complete models could be an interesting topic of further investigation, and could allow our 5d simplified models to be used in full event generators for experimental studies. In this regard, it would be interesting to determine the effect on event shape variables of replacing our massless particles, which stand in for SM particles in the current study, with QCD jets from quarks and gluons, or with relatively soft photons. In the LHC context, one would need to consider carefully the impact of initial state radiation, the underlying event and pileup. It is also not clear what event-shape variables would be most effective in reducing backgrounds at a hadronic collider. All of these issues must be addressed before optimized searches for phenomena of this type can be designed.
In particular, notice that if y = 0, only terms with n = 0 contribute, and the formula becomes independent of d and reduces to an ordinary hypergeometric function 2 F 1 : Similarly, F 4 (a, b; c, d; 0, y) = 2 F 1 (a, b; d; y).
There is a general result that relates an infinite integral of products of three Bessel functions times a power to the function F 4 [86]: valid (at least) when the λ, ν and m parameters are real, m 1 > m 2 + m 3 , λ + i ν i > 0, and λ < 5/2. Notice that the couplings we wish to calculate, (4.1), have the form of this integral in the special case λ = 2. On the other hand, identifying the first four arguments of the F 4 function as α, β; γ, δ, we see that they obey α + β = γ + δ + λ − 2. (A.7) Thus, the reduction (A.5) to ordinary hypergeometric functions applies when λ = 1, which is not the case of our interest. However, integrals with different values of λ may be related to each other using the identity This allows us to obtain the λ = 2 case of our interest from the λ = 1 case with a known reduction by taking a derivative with respect to a mass parameter. This method of relating different integrals was discussed in [94,95], and a specific formula relevant for the case of our interest was given in [96]. This approach leads to the equation (4.13) that we have given in the main text.

B Dependence on n p
In this section we explore the event shape dependence on the initial mode n p with Dirichlet boundary conditions. As we increase n p , this can affect both the coupling structure of the initial decays (where the degree of KK-number violation affects the event shape) and the kinematics of the final state.
To understand how the kinematics of the final state depends on n p , we first reconsider the spherical toy model described in section 3.1, i.e., events with n p identical particles at rest, which split into n p pairs of massless particles with equal energy and random orientation. For sufficiently large n p , the randomly distributed particles will appear isotropic. In figure 17, we see that as n p increases, the distribution of event isotropy I sph 192 converges toward the theoretical limit of approximately 0.1, by eq. (1.5).
Similar behavior is seen in figure 17 for the single field example with ν = 0, which has slightly larger event isotropy than the spherical toy model due to non-zero KK-number violation and energy anisotropy. More generally, we expect near-spherical examples to converge to the theoretical bound at large n p .
To explore the highly non-spherical regime, we consider another toy model that produces boosted particles early in the cascade. We generate a back-to-back pair of initial particles of equal mass m = np 2 m 0 , where m 0 is the mass of the lightest HSH, and with boost γ. Each particle decays to n p /2 HSHs which are at rest in the particle's decay frame. The HSHs then decay to two massless particles each. Thus the two initial particles produce two quasi-spheres, made of n p massless particles with equal energy, that are boosted by γ in opposite directions.
We study the event shape dependence on n p as we increase the boost γ. When the initial particles are produced at rest (γ = 1), we recover the spherical toy model, as in figure 17, where the isotropy gradually decreases with n p . In the opposite limit of large If I + vanishes, which occurs when ν 2 + ν 3 − ν 1 is a positive even integer, then there is a subtlety. The partial width of a typical decay dominated by I − ∼ √ m 1 m 2 m 3 λ −3/2 PS scales as n −3 p . This would suggest that IT region dominates and that the total width scales as n −1/2 p . While the latter conclusion is correct, near-threshold decays, with small kinematic boosts but larger KK-number violation, can be as important as the minimally-KK-violating decays at the IT. There is therefore a band in the branching fraction triangle, including but extending beyond the immediate threshold, which remains important as n p → ∞. The boosted region scales away, so the events are far from jetty. However, because KK-number violation is non-minimal, the events may not become spherical in the n p → ∞ limit.