Unsupervised Hadronic SUEP at the LHC

Confining dark sectors with pseudo-conformal dynamics produce SUEP, or Soft Unclustered Energy Patterns, at colliders: isotropic dark hadrons with soft and democratic energies. We target the experimental nightmare scenario, SUEPs in exotic Higgs decays, where all dark hadrons decay promptly to SM hadrons. First, we identify three promising observables, the charged particle multiplicity, the event ring isotropy, and the matrix of geometric distances between charged tracks. Their patterns can be exploited through a cut-and-count search, supervised machine learning, or an unsupervised autoencoder. We find that the HL-LHC will probe exotic Higgs branching ratios at the per-cent level, even without a detailed knowledge of the signal features. Our techniques can be applied to other SUEP searches, especially the unsupervised strategy, which is independent of overly specific model assumptions and the corresponding precision simulations.


Introduction
Hidden sectors are one of the most interesting and generic paradigms for physics beyond the Standard Model (BSM). These kinds of new particles and forces are not only plausible from a bottom-up point of view, but also arise in many top-down BSM theories . In most scenarios of interest, hidden sector particles couple to SM fields via feeble interactions or heavy messengers, and the nature of these portal couplings determines their collider phenomenology. A case of special interest are hidden valleys [1], referring to hidden sectors with a confining gauge group which gives rise to rich infrared (IR) dynamics from very simple ultraviolet (UV) theory structures. The production of dark quarks at the LHC leads to a dark shower and high-multiplicity production of dark hadrons, in analogy to QCD jets. Depending on the portal by which the dark hadrons are produced and decay, these dark showers produce a wide variety of LHC signatures, which have been the subject of intense theoretical and experimental study in the last decade [2,8,[25][26][27][28][29][30][31][32][33][34][35].
We consider one of the most challenging varieties of dark showers, Soft Unclustered Energy Patterns (SUEPs). If a hidden valley possesses a large gauge coupling that is pseudo-conformal above its confinement scale, then large-angle emission is unsuppressed for most of the parton shower evolution. This means the dark hadrons are not arranged in narrow QCD-like jets, but emitted approximately isotropically in the shower centre-ofmass frame [36,37]. This defines the SUEP final state as a high-multiplicity sphericallysymmetric shower of hidden sector states. While existing searches can be sensitive to SUEPs produced at high energy scales, or with dark hadrons that decay to sufficiently conspicuous final states like leptons or Long-Lived Particles (LLPs) [8,37], no dedicated SUEP searches exist to date. Furthermore, SUEPs without conspicuous final states are not captured by any existing search and represent an unusually cruel signal, since their soft, isotropic distributions can mimic the ubiquitous pile-up produced by simultaneous LHC collisions.
To ensure that all types of SUEP signals can be discovered at the LHC, we focus on a well-motivated SUEP nightmare scenario, where SUEP is produced in exotic Higgs decays and the dark hadrons decay promptly and exclusively to SM hadrons. The modest energy scale of exotic Higgs decays and the lack of conspicuous final states forces us to rely on the kinematics of the resulting SM hadrons to extract the signal from the overwhelming QCD background. This production mode also allows us to side-step the problem of how to trigger on SUEPs [37] by using leptons in associated V h production. The analysis techniques we develop will not only allow for the detection of this SUEP nightmare scenario, but should also increase the LHC experiments' sensitivity to all other SUEP possibilities.
An acute obstacle for SUEP searches is the lack of rigorous predictions and simulations for the strongly coupled pseudo-conformal dark sectors. Rather than hoping for conspicuous final states, we utilize the kinematics of the SM hadrons without relying on fine details of the signal beyond the robust SUEP characteristics of isotropic, soft, and democratic dark hadron energies. We therefore simulate SUEP production using a simple QCD fireball model of thermal dark hadron emission [20], and find robust observables and event representations. An important observable for SUEP searches is the inter-particle distance matrix ∆R ij between charged hadrons. This matrix encodes the essential geometric differences between QCD-like and SUEP-like hadron production. It forms the backbone of all our analysis strategies, together with known variables like event isotropy [38] and the charged particle multiplicity.
To demonstrate the distinguishing power of our observables, as well as the drastic improvements from more sophisticated techniques, we examine three strategies for systematic SUEP searches at the HL-LHC. First, we simulate a simple cut-and-count analysis, which will turn out to allow for impressive sensitivities to Br(h → SUEP) ∼ 1% at the HL-LHC. We anticipate that a realistic analysis with data-driven background estimation will perform even better, since our study is limited by background simulation statistics. To improve on this, we utilize supervised as well as unsupervised machine learning (ML). Unsupervised analysis concepts along the lines of autoencoder neural networks [39] have the potential to transform LHC analyses [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56], including searches for dark showers [45]. We point out how unsupervised methods are especially appealing for difficult-to-simulate signals like SUEP since they only rely on the known QCD background for training, while yielding at least several times greater SUEP sensitivity than the cut-and-count analysis.
Our investigation establishes that SUEP searches need not rely on conspicuous SM final states for excellent sensitivity at the HL-LHC. The unique dark hadron kinematics, which robustly follows from their origin in pseudo-conformal strongly coupled dynamics, allows for the SUEP final state to be distinguished from its overwhelming QCD background. Our techniques can be applied to all SUEP searches to dramatically enhance their sensitivity, regardless of energy scale or SM final state. This paper is structured as follows. In Section 2, we briefly review the SUEP theory and define the benchmark scenario for our study. Signal and background simulation is discussed in in Section 3. In Section 4 we define the relevant observables to distinguish SUEP from QCD and discuss supervised and unsupervised machine learning techniques. In Section 5 we present our results, including projections for the Br(h → SUEP) sensitivity at the HL-LHC, and we conclude in Section 6.

Theory of SUEPs
Hidden valley models are a large class of BSM theories in which the SM is extended by additional gauge groups under which SM particles are neutral. New particles that are charged under non-Abelian extensions can give rise to a wide range of interesting hidden sector dynamics [1,2,25] and various challenging SM signatures. These models appear in the context of many top-down constructions [57], including string theory, of course [21], and are compatible with various potential resolutions to the hierachy problem such as supersymmetry [58], little Higgs models, TeV extra dimensions and Randall-Sundrum scenarios [1,[59][60][61]. In Neutral Naturalness scenarios, hidden valleys actually solve the little hierarchy problem directly [3]. Models containing a hidden valley have also been studied in the context of dark matter [22], matter-antimatter asymmetry [23] and the origin of neutrino masses [24].
In most hidden sector scenarios with a confining gauge force, the dark parton shower is qualitatively similar to QCD: the asymptotic freedom of the running gauge coupling enhances soft and collinear emission, resulting in the production of hidden sector states in collimated jets. If the hidden 't Hooft coupling is large (λ ≡ g 2 N c 1) and approximately constant over a significant energy range, the distribution of the produced dark hadron momenta will be much more democratic than the hierarchical jet-like behaviour we see in QCD. This is the SUEP class of signals, characterized by relatively soft, isotropic emission of dark hadrons [36,37]. It is worth keeping in mind that strongly coupled hidden sector dynamics is not the only scenario that can lead to SUEPs. Similar final states can be produced in many-step cascade decays in the hidden sector [62], or in theories with extra spatial dimensions (see e.g. [63][64][65]). We now describe the general features of production, evolution and decay that constitute the SUEP signal, and define a benchmark scenario for our study.

Dark Hadron Production in Exotic Higgs Decays
Production of SUEP can occur through various portals coupling SM particles to SM-singlet states charged under a dark gauge group. Higgs and vector boson portals are commonly studied examples. Alternatively, a new particle charged under both SM and the dark gauge group could be produced [20]. We assume that hidden sector states ψ D can be produced via the Higgs portal: O production ∼ |H| 2 ψ DψD .
(2.1) ψ D could represent a fermion or scalar dark quark charged under the hidden gauge group.
In the fermion case, the above operator is actually dimension 5 with a coupling ∼1/Lambda for some UV scale Λ. This operator will give rise to the production of hidden states in exotic Higgs decays, as long as the dark hadrons are kinematically accessible. The dark quarks hadronize into a large multiplicity of dark hadrons. For a confining gauge theory, the evolution of the system from the hard scale Q at which the first partons charged under the gauge group are generated to the IR confinement scale Λ governs the hadron multiplicity generated during showering. Q is of the order m h in our case. For QCD jets, this evolution can be reliably simulated as a parton shower, but for SUEPs different strategies are required.
The average hadron multiplicity is given by where F is the fragmentation function describing the final state distribution of momenta [66], and x is the momentum fraction of a given splitting. Calculating its evolution from the scale Q down to an IR scale Λ involves resumming the divergent contributions from the anomalous dimensions, with the leading contribution obtained from just the first (j = 1) Mellin transform of the anomalous dimension. If the theory is conformal between Q and Λ, the running of the coupling can (by definition) be neglected. In the limit where the 't Hooft coupling is large λ (1 − x) it was first shown in Ref. [67] that where γ T (1) is the first Mellin-Moment of the time-like anomalous dimension of the fragmentation function. Further expanding to zeroth order in the coupling yields Note that the small momentum fraction x carried by each individual splitting follows from λ (1 − x). In this strong regime branching is expected to yield emissions with x ∼ Λ/Q that are relatively isotropic in direction and democratic in momentum [37,67]. Thus, with a large enough scale separation Q Λ, low-x, high-multiplicity final states are generated. Branching terminates after N final ∼ log n splittings at Q N final ∼ Q/2 N final ∼ Λ, at which point hadronization takes over.
Hadron production in QCD at high temperatures is a close analogue to our situation, and statistical models have consistently shown that hadron multiplicities follow a thermal distribution [67][68][69][70][71][72][73][74][75][76][77]. We use this picture as a toy model of dark hadron production in SUEP, modelling the distribution of dark meson momenta as a relativistic Boltzmann distribution where m D is the mass of the final dark states and T D acts as the Hagedorn temperature of the hidden confining gauge force, with T D ∼ Λ [37,71]. This temperature controls the kinetic energy of the dark hadron distribution.

Dark Hadron Decay
Generically the decay of the dark hadrons φ D into the SM will occur through some effective coupling of the form where O SM contains fields charged under the SM gauge group. The phenomenology of the SUEP signal is determined to a large extent by the dominant decay portal. For example, if the dark hadrons decay to massive dark photons that mix with the SM photon, then the SM final state contains both hadrons and leptons in roughly gauge-ordered proportions, though for various dark photon masses in the GeV-regime, hadrons can dominate [37]. Alternatively, if the dark hadrons decay through the Higgs portal, the final state will contain hadrons and leptons in roughly Yukawa-ordered proportions. For sufficiently small portal coupling, the decay length of the dark hadrons may also be macroscopic, resulting in LLP signatures [8].
It is also possible for the hidden sector states to decay purely hadronically, which is the most experimentally challenging case at the LHC. A very simple example is the gluon portal as described in [20]: Here, a is a heavy elementary pseudo-scalar in the dark sector and ψ D is the dark quark. Dark hadrons, which are bound states of ψ D , could then decay to SM hadrons via an effective operator φ D GG. Another example is the hadrophilic (or leptophobic) Z portal [78], where a new heavy gauge boson couples to SM quarks but not leptons, allowing dark hadrons to decay via an effective operator like φ D qq.

Prompt Hadronic Benchmark Scenario
The production of hidden sector states via the Higgs portal generally and exotic Higgs decays in particular is one of the most motivated and plausible discovery scenarios for new physics [79]. It is therefore vital that our experimental search strategies cover all possibilities for a signal at the LHC. This is especially urgent since for final states that are not covered by existing searches, branching fractions of ∼10% are easily allowed by current measurements of Higgs couplings and invisible decays [80][81][82][83][84].  Figure 1: Cartoon of the our benchmark scenario for SUEP produced in Higgs decays with prompt decay of dark hadrons into SM hadrons. The SUEP description applies in the purple area, for T within a factor of a few of the dark hadron mass. In the green region, dark hadron production is not thermal, but described by processes more analogous to chiral QCD. In the orange and green regions, searches for Higgs decays to a few resonances would be sensitive to this dark sector. The parameter space region indicated with the darker purple rectangle is the focus of our analysis. Our cuts are optimized for the most archetypically SUEP-like final states, schematically indicated by the lower-left corner of this rectangle, demarcated with the diagonal line.
We therefore focus on the experimental worst-case scenario for SUEP produced in exotic Higgs decays: purely hadronic and prompt decays, with a particular interest in low dark hadron masses that make resonance searches [85,86] or applications of jet substructure techniques [87] challenging. While the simplest gluon portal scenarios suggest that dark hadrons lighter than ∼10 GeV have macroscopic decay lengths [20] (which could allow for the use of long-lived particle search techniques [8,33]), other possibilities can easily realize prompt, purely hadronic decays over a much wider range of dark hadron masses. For example, the hadrophilic (leptophobic) vector portal [78] with a hypothetical confining sector where the lightest dark hadron is a dark-rho-like vector ρ D would allow dark hadrons lighter than a GeV to decay promptly into SM hadrons. In focusing on the prompt case we can develop techniques that allow SUEP production to be identified using only the geometrical and momentum distribution of its SM final states. These techniques will enhance our sensitivity for any kind of general SUEP, in addition to whatever other features of the final state, like displaced decays, leptons, or photons, can also be exploited.
The parameter space of our benchmark scenario is shown schematically in Figure 1. The Higgs production portal sets the high scale for the event to Q = 125 GeV. With this scale fixed, the distribution of final states is in principle determined by the dark hadron mass m D and the dark Hagedorn temperature T D of Eq. (2.5), shown on the vertical and horizontal axes of Figure 1. In reality, there may be different dark hadrons with different spins, decays, and distributions, but this simplified description is sufficient for our purposes. On the left and right, the relevant parameter space is kinematically bounded by the red-hatched areas. Dark hadron production in Higgs decays requires m D < m h /2. As we explain below, we focus on dark hadrons that decay to SM hadrons, which requires m D > (2 − 3)m π , depending on the exact decay portal.
Even within that mass range, hidden sectors with pseudo-conformal dynamics do not always manifest as SUEP signatures in exotic Higgs decays. If m D > m h /3 (blue area in Figure 1), then dark hadrons are only pair-produced in Higgs decays, making this scenario equivalent to standard exotic Higgs decays to pairs of various new particles (see e.g. [79]). If the dark hadron mass is fairly large, m D ∼ m h /(few) (yellow area), or the dark Hagedorn temperature is comparable to or above the Higgs mass T m h (green area), then exotic Higgs decay would produce only a small multiplicity of dark hadrons that are either fairly hard or fairly heavy. In both cases, dark hadron production is not thermal but is described by processes more akin to chiral QCD. These regions are likely accessible through modified searches for resonances in exotic Higgs decays [85,86], and we do not focus on them here. The gray area where T /m 1 is not expected to be realized by any pseudo-conformal hidden sector, since the dark hadron mass and temperature are both related to the strong coupling scale Λ. On the other hand, T /m 1 is possible in what we call the "dark pion regime", where the dark hadrons are pseudo-Goldstone bosons of an approximate symmetry, meaning their mass can be much smaller than Λ. We do not focus on this region, but it would be an interesting target for future investigations.
This leaves us with the actual SUEP regime for dark hadron production in exotic Higgs decays, indicated by the light purple area. In this investigation, we will focus on dark hadron masses below 8 GeV and T /m in the reasonable range of ∼ 0.25 to 4. This target SUEP parameter space is marked out as the darker purple rectangle. Our cuts will be particularly optimized for the lower-left region of the rectangle, demarcated with the thick diagonal line. This is the region of low dark hadron mass and/or temperature, corresponding to the softest SM final states that are most difficult to search for using existing techniques.

Simulation
We briefly outline how we simulate our SUEP signal and the most important QCD backgrounds, where the latter is necessary to develop our analysis techniques, even though a realistic experimental analysis would rely on data-driven background estimation.

Signal
We generate event samples for exotic Higgs decay into SUEP using the SUEP Generator plugin [20] for Pythia 8.243 [88], which models the dark shower as a spherical distribution of dark pseudo-scalar mesons with momenta drawn from the relativistic Boltzmann distribution Eq. (2.5). As in Ref. [37], we make the simplifying assumption that there is only one flavor of dark meson produced in the exotic Higgs decay. The free parameters of the SUEP shower are the dark hadron mass m D and the effective temperature T D , setting the energy scale at which dark hadrons are produced dominantly.
We simulate associated Higgs production at the 14 TeV HL-LHC, pp → V h, V → / ν, in Pythia 8. The Higgs is decayed to the SUEP final state of dark mesons, which then decay directly to a uū quark pair that in turn undergoes SM hadronization. The exact choice of hadronic decay mode does not significantly affect our analysis, so we use this single channel as a stand-in for other purely hadronic portals. The events are then passed through the simplified detector simulation code Delphes 3 [89] with CMS detector settings. The simulated detector-level objects output by Delphes are used for our analysis. Our SUEP search will only use charged-track information. Since charged tracks can be traced back to the primary vertex, they are very robust with respect to pile-up contamination. We can therefore neglect the effects of pile-up in the remainder of our study.
In our event samples we cover m D from 400 MeV to 8 GeV, and T D /m D from 0.25 to 4. The lower bound on the dark hadron mass ensures that the dark mesons are kinematically allowed to decay to two pions. The upper bound is chosen to show where our search loses sensitivity. The range of temperatures is chosen to satisfy T D ∼ m D , which is the regime where the thermal picture of SUEP production is valid. The signal cross section is [90] for SUEP production in exotic Higgs decays in association with leptonic W/Z-bosons that decay into electrons or muons. In total, we generate 2.0×10 5 Zh and W h signal events, proportional to their respective cross section, for each set of signal parameters (m D , T D /m D ).

Background
The dominant background to our signal is production of one or two leptons in association with any number of QCD jets. It is highly challenging to model reliably, and in a realistic study, data-driven background estimation would be employed, see Section 4.4. However, for the purpose of developing our analysis techniques, we simulate QCD+leptons background samples using MadGraph5 aMC@NLO 2.6.6 and Pythia 8.243. Ideally, one should simulate fully matched multi-jet + / ν samples to capture the background distribution as closely as possible. However, due to the large statistics needed for our analysis, and the fact that such a simulation is anyway unlikely to be a perfect representation of the detailed hadronic distributions at the relevant high multiplicities and relatively low energy scale of the Higgs mass, this approach is not practical. Instead, we simulate nj + / ν, where n = 2, 3, 4 without jet matching and p T > 15 GeV at generator level, to determine the effect of jet multiplicity at the hard event level on our analysis. We find that n > 2 leads to lower cross section while being more distinguishable from the SUEP signal using the analysis techniques we develop here. Therefore, to be conservative, we simulate 2j + / ν as our background samples for Zh and W h production and decay into SUEP, respectively. In total, we use 10 8 background events to represent the lowest-order MadGraph5 cross section for this background. While this is sufficient to develop our analysis techniques, the Monte Carlo background sample has ∼ 1/100 the statistics of the full HL-LHC dataset. This is important for the interpretation of our results in Section 5.

Analysis
The goal of this paper is to devise strategies for extracting SUEP signals from a large background events without relying on the details of the simulated signal. We first describe our trigger assumptions and baseline cuts, and define a cut-based classifier to establish how sensitive such a simple approach can be. Section 4.2 introduces a supervised neural network classifier, to demonstrate both the advantages and limitations of the supervised approach for our physics problem. We then introduce our primary tool in Section 4.3an unsupervised neural network that we employ as an anomaly detector for SUEP.

SUEP Observables
We define our trigger pre-selection by requiring that all events have at least one charged electron or muon with p T ≥ 40 GeV, or two opposite sign charged leptons with p T ≥ 30 (20) GeV. We also require that the scalar p T -sum of hadronic charged tracks from Delphes is above 30 GeV. Both signal and background have a trigger efficiency of ≈40%, relative to the cross sections in Eq. (3.1) and Eq. (3.2). We focus on the most challenging region of SUEP parameter space, with either low dark hadron masses or low dark shower temperatures. This gives rise to the most archetypically SUEP-like final states with a high multiplicity of isotropically distributed, relatively soft SM hadrons. The three observables that best capture the characteristics of this signal are the charged particle multiplicity N charged , the event isotropy I [38], and the interparticle distance. In all steps of the analysis that follow, we only use charged particle tracks with p T ≥ 300 MeV from the primary vertex, excluding the one or two hard leptons associated with the decaying gauge boson.
The event isotropy observable I ∈ (0, 1) [38] quantifies the energy mover's distance between a collider event and an idealized isotropic event with uniform energy distribution, so I = 0 indicates a fully isotropic event. Originally, three different definitions of the event isotropy are laid out, utilizing different geometries -spherical, cylindrical, or in a two-dimensional ring. We compute the ring isotropy of the set of charged hadronic tracks of each event, since at a pp collider we have no way to know the longitudinal boost of the Higgs that decays to SUEP. Since the SUEP is isotropic in the Higgs rest frame, we boost the hadronic charged track system of each event into its transverse rest-frame before computing the ring isotropy. To do this we assume that all hadronic charged tracks belong to particles with the pion mass, but this is sufficient to significantly separate signal and background events. The variable that we introduce for the specific purpose of studying SUEP events is the interparticle distance matrix ∆R ij for charged hadron tracks in the lab frame. It captures the unique topology of SUEP events while being very suitable for machine-learning applications, since it is invariant under re-definitions of the azimuthal angle around the beam axis. It is also useful to define the mean ∆R of all matrix entries. Figure 2 shows distributions of N charged , I, and ∆R for the QCD background and a variety of SUEP benchmark points after trigger selection. The separation between signal and background is clear, with SUEP having higher multiplicity, more isotropic distribution of tracks, and a significantly wider spread of inter-particle distances. To understand how these observables change across SUEP parameter space, we show their average values (across the whole sample after trigger selection) as a function of m D and T D in Figure 3. The pairwise correlations between each of these observables are included in the Appendix in Figure 8, demonstrating that each of these three variables encode distinct information about each event.
In addition to the trigger requirements, we therefore impose the following pre-selection cuts on all events: N charged ≥ 70 , I < 0.07 , ∆R < 3 . This cut targets the most SUEP-like parts of signal parameter space with low dark Hagedorn temperature and/or low dark hadron mass (see Fig. 1). All but 2.2% of the posttrigger background is eliminated, while leaving 31.8% of signal for m D = 0.4 GeV and T D = 0.4 GeV. These requirements are less optimal for larger dark hadron masses or temperatures -for example, the signal efficiency is only 1.1% for m D = 5 GeV, T D = 20 GeV. However, larger temperatures and masses generally lead to higher-energy final states or separable resonances, and are therefore not the focus of our present analysis.
We now consider three options for SUEP searches at the HL-LHC, using events which pass the baseline pre-selection as a starting point: 1. The simplest strategy is a cut-and-count analysis using high-level observables. It will serve as a baseline for more sophisticated machine-learning techniques, and can be implemented very easily and effectively with a stricter cut on ∆R compared to Eq. (4.1). Varying the cut threshold yields a significance improvement curve (SIC) of the signal efficiency vs the background efficiency for each point in signal parameter space, which we will then be able to compare to results using machine learning methods. As we will see, this already yields very promising SUEP sensitivities.

2.
A supervised ML-classifier requires detailed knowledge of the signal, since it is trained on signal and background. This makes supervised techniques unlikely to be a realistic analysis candidate for broad SUEP searches. However, we perform a simple supervised study in Section 4.2 to demonstrate the best-case scenario for SUEP sensitivity if the signal was very well understood.
3. In Section 4.3, we use an unsupervised autoencoder trained only on the background as an anomaly detector to improve on the sensitivity of the cut-and-count analysis. This is likely to be a realistic analysis candidate since it can be performed using datadriven background estimation techniques without precise knowledge of the signal.

Supervised ML-Classifier
While a supervised classifier is explicitly dependent on the characteristics of the signal (and background) simulation used in its training, which we cannot trust in detail, it can still provide a useful comparison to an unsupervised network, and give an indication of how sensitive a search based on similar methods could be with improved modelling of the SUEP signal. An additional limitation is that even with reliable signal simulation, the parameters of the real SUEP are unknown, and a supervised network trained to recognize SUEP with one set of parameters may fail if the true parameters change. As we will see, this is indeed the case. The supervised network architecture we choose is a dynamical graph convolutional neural network [91,92]. We implement the network in PyTorch. The input feature representation for both the supervised and unsupervised networks is the interparticle distance matrix, with the redundant left-lower half set to zero and track p T information added to the diagonal: All events are required to have at least N charged ≥ 70 tracks, and all events are truncated to keep only the 70 highest-p T tracks to ensure each event has the same dimensionality in the analysis. The input matrix ∆R ij is used to generate graph edges between each particle (node) and its k = 7 nearest neighbours in ∆R space. The node features for each particle are the 70-dimensional vector of ∆R distances to all other particles in the event. The graph network has two EdgeConv blocks, each comprising a three-layer perceptron with leaky ReLU activation [93]. The EdgeConv operation updates the node features x i of each particle as where k is the number of neighbours assigned to each node, and h θ is a non-linear function of learned parameters θ, implemented as a three-layer perceptron. The graph edges are re-computed between the first and second blocks using the Euclidean distance between the feature vectors of each node to determine its 14 nearest neighbours. The first block has feature dimension 64, while the second block has feature dimension 32. Batch normalization follows each layer. The output of the graph layers is averaged, then passed through two fully connected layers, first expanding to dimension 128, then down to output dimension 2. The loss function is the cross-entropy loss between the output and the true class label of each event.
L(x, class) = − log To test how the performance of the model depends on the choice of training parameters, we train twelve neural networks on twelve different choices of (m D , T D /m D ), with m D = 0.5, 1, 2 GeV and T D /m D = 0.5, 1, 2, 4, and evaluate their effectiveness over the whole signal parameter space. We also evaluate the efficacy of a 'cocktail approach' [92,[94][95][96] by training on a mixed sample including signal events from each of the twelve parameter choices. A conditional training [97,98] on the signal parameters would be a similar approach. Each network is trained for 10 epochs with a decaying learning rate. Longer training periods were tested and found to be unnecessary. Loss values for each test sample event are obtained for the model realizations of the last 5 training epochs before being averaged.

Unsupervised Autoencoder
Supervised ML-classification is an extremely powerful analysis tool, but it only works for signals with well-defined and universal features and corresponding precision simulations. At the LHC, this is not always the case, and SUEP with its toy shower is a perfect example for a more broadly defined signal. Here we prefer not to train a classifier on signal simulations. Instead, we employ anomaly detection methods, where we train a network only on the well-understood background dataset, so that it can flag events that are anomalous in comparison. We use an autoencoder [39,45,46] and train it on data without class labels, with a loss function that incentivizes its output to be as close as possible to the input. The intermediate network layers have a restricted number of nodes compared to the dimension of the input and output, forcing the network to compress the information in the input, and then decompress it to recover the output. The principle of the autoencoder's use as an anomaly detector is that it should fail to accurately reconstruct events that are anomalous compared to the dataset it was trained on. A high reconstruction loss flags an event as being potentially anomalous, or in collider physics parlance, a signal (SUEP) candidate event.
The autoencoder uses the same modified ∆R ij matrix input as used by the supervised network, see Eq. (4.2). Other representations were tested, including the matrices of both k T and anti-k T distances d ij = min(k ±2 ti , k ±2 tj )(∆R ij ) 2 /R between particles [99], the high-level observables used in the pre-selection cuts, as well as the raw p T , φ, η values for each particle; all yielded results that were much less useful than the analysis we present here. This emphasizes the importance of choosing the correct input representation over sophisticated network architecture for SUEP searches.
The matrix ∆R ij is flattened into a vector of length N 2 charged and fed into the autoencoder. The neural network comprises five fully connected layers, with the number of nodes decreasing to the bottleneck size in the third layer, then increasing back to N 2 charged = 4900. An alternative 3-layer network performs only slightly worse. For the bottleneck size we find that larger bottlenecks consistently lead to better performance than smaller ones, so we use N bottleneck = 1000. Each layer of the network has a leaky ReLU activation with slope -0.2 for negative values of x, except the final layer which has a ReLU activation. The loss function of the network measures the difference between the network's input and output as where σ(x) = 1/(1 + e −x ). Among different values of m and n, starting with the usual mean squared error m = 2, n = 0, the best-performing is m = 3, n = 0. The sigmoid normalization of the input is essential to the network's success. Without it, the autoencoder encodes SUEP events with slightly lower loss than the QCD background on which it was trained, and completely fails to identify anomalous events. We hypothesize that this is because, unlike many experimental signatures, SUEP is less complex than its QCD background, and it has smaller values of ∆R than QCD. This complexity bias has been noted before [45,52,53]. The sigmoid function reduces sensitivity to large values of ∆R and p T by mapping them to values very close to 1, while remaining approximately linear for small input values, but offset to a minimum value of 0.5. These effects make it easier to accurately reconstruct QCD events while enhancing the network's sensitivity to deviations from the input on the SUEP events with characteristically smaller absolute values of the input features. Out of 8.8 × 10 5 background Monte Carlo events that pass the pre-selection cuts, 2.4×10 5 are used for training, 5×10 4 for validation when tuning network hyperparameters, and 5.9 × 10 5 for testing. The number of signal events that pass the cuts varies with m D and T D , but generally a few ×10 4 events remain at each parameter point to be used for testing. SUEP events with m D = 1 GeV, T D = 0.5 GeV are used for validation purposes.
The network is trained for 15 epochs with a decaying learning rate. Longer training periods were tested and found to be unnecessary. Loss values for each test sample event are obtained for the model realizations of the last 5 training epochs before being averaged.
Other architectures were investigated, including variational autoencoders [100] and a graph convolutional autoencoder utilizing the same EdgeConv operations as the supervised network, which itself led to the use of the ∆R ij event representation. Interestingly, the simpler, fully connected architecture consistently delivered much better background rejection than any of the more sophisticated graph networks in the unsupervised approach.

Data-driven Background Estimation
Following the logic of this section further, we briefly describe how a data-driven background sample could be derived for use in a realistic experimental analysis based on our study. The total background cross-section can be measured directly (and compared against simulation), since the leptons+QCD signal region is completely background dominated for SUEP production in exotic Higgs decays. The background efficiency of the classifier, whether it is based on cuts, an unsupervised neural network, or a supervised one, can be estimated by evaluating it on a control region.
A variety of control regions are possible, but perhaps the most promising is defined by replacing the lepton criterion by a mono-photon criterion. Such a sample should be free of signal contamination, and its hadronic content is extremely similar to the hadronic background to our SUEP search, since it is unconstrained in its detailed production channel except that it recoils off a hard electroweak boson. The QCD distributions in this monophoton-plus-jets sample should therefore closely match those in the Z+jets signal region, especially if the control region sample is reweighted to match the shape of the photon p T spectrum to the dilepton p T spectrum in the Z+jets signal sample. A variant of this approach can likely be adapted to estimate the background in the W +jets channel as well, by using simulation to compute a transfer function from p T,W to p T, and applying it to p T,γ in the control region before applying the p T, reweighing.
Based on the assumption that such a strategy can be implemented, we therefore will use benchmark estimates of 1% and 10% for the systematic background uncertainty in our analysis to estimate the final physics reach.  . For the autoencoder, the peak sensitivity improvement we are able to probe reliably with the statistics of our Monte Carlo samples (see text) corresponds to signal and background efficiencies of (1.6×10 −2 , 1×10 −7 ), (3.7×10 −4 , 1.5×10 −7 ), (6.1×10 −2 , 1×10 −7 ) and (2.2×10 −2 , 1×10 −7 ) relative to trigger selection, respectively. efficiency (after triggering) or number of remaining signal events at the HL-LHC as the threshold is varied. The vertical axis shows signal significance improvement relative to the trigger sample. From these specific examples, a few general features are apparent:

Results
• The autoencoder consistently outperforms the simple ∆R cut significantly; • The supervised networks outperform the autoencoder for signal parameters close to their training values, but can perform much worse for different parameters; • Our analysis is not optimized for larger dark hadron masses and temperatures; • For lower dark hadron masses or temperatures, both the cut-and-count and autoencoder analysis strategies are very powerful, yielding orders of magnitude improvement in signal significance compared to the baseline preselection cuts of Eq. (4.1).
In fact, the statistics of our QCD background sample is insufficient to capture the true power of our analysis techniques. We would expect the significance improvement to increase as the signal efficiency decreases, but only until the SI curve turns around when the cut becomes too harsh. This is clearly the case in the top-right panel of Fig. 4, for signal parameters to which our analysis is not optimized. However, for the three other examples, we never reach this turn-over point before running out of simulated background events. It is not even clear if this turn-over is reached with the full statistics of the HL-LHC (100 times greater than our simulated background sample).
To understand how much harder we might be able to cut on the background, Fig. 5 shows the same SIC curves but with background efficiency on the vertical axis. Naively extrapolating, we can anticipate that a realistic autoencoder search with a fully data-driven background sample at the HL-LHC might be able to reach the very-low background regime while retaining enough signal to probe Br(h → SUEP) as much as an order of magnitude smaller than the sensitivity estimates we can rigorously derive here.
As a result, the reach projections we present in this paper will be very conservative. Furthermore, the limited statistics of our background sample means that reach projections will be very similar for the three analysis methods despite their obvious differences in performance in Figs. 4 and 5. It is therefore important to additionally evaluate our classifiers using a somewhat orthogonal metric.
The area under the curve (AUC) of the signal efficiency as a function of background rejection can be computed as a function of the number of events kept after the baseline preselection cuts, see Eq. (4.1). This metric is standard to measure the performance classifiers independently of the choice of threshold. Fig. 6 shows the AUC achieved by the cut-based classifier, the fully connected autoencoder, and the supervised graph network trained on a cocktail of signal events, where for the latter the signal parameters are indicated with red dots in the (m D , T D /m D )-plane). Unsurprisingly, the highest AUC values are achieved at low m D or T D /m D , with the supervised graph network modestly outperforming the autoencoder, which outperforms the simple cut, across the SUEP parameter space. The performance of supervised networks trained on single signal parameter points are presented in the Appendix.
A key physics result is the actual sensitivity of HL-LHC searches to the SUEP final state. We therefore extract the smallest Br(h → SUEP) branching ratio for which S/ B + u sys B 2 > 2, where u sys gives the systematic uncertainty on the background. Because of the high degree of background rejection required to be sensitive to relevant Higgs branching ratios, and the limited size of our background dataset, very few or no simulated background events remain when the classifier's threshold is set to maximize sensitivity to low branching ratios. This introduces a significant statistical uncertainty to the estimated LHC reach. We can decouple our estimate somewhat from these limitations by demanding the statistical uncertainty from the limited size of the background dataset to remain below 50%. This is conservative, since a realistic analysis will employ even harsher cuts with better significance improvement. For the same reason, the difference in reach between our three methods is likely underestimated. The performance gaps between the ∆R cut, autoencoder, and supervised network are more clearly shown by the individual significance improvement curves and the AUC differences.
Finally, Fig. 7 shows the SUEP sensitivity achievable by the ∆R cut, the autoencoder, (b) (c) Figure 6: AUC for (a) cut on ∆R, (b) fully connected autoencoder, (c) supervised graph network trained using cocktail of signal parameter choices (training parameters indicated with red dots). These plots illustrate the significant performance improvements of the autoencoder relative to the simple ∆R cut, and of the supervised cocktail approach relative to the autoencoder.
(a) (b) (c) Figure 7: Minimum excludable Br(h → SUEP) at the HL-LHC, assuming 1% systematic uncertainty on QCD background for (a) cut on ∆R, (b) fully connected autoencoder, (c) supervised graph network trained using cocktail of signal parameter choices (indicated with red dots). Note that the limited statistics of our QCD background sample leads these projections to be very conservative while also de-emphasizing the performance differences between the three methods. and the cocktail-trained supervised network, all assuming 1% systematic uncertainty on the background. Both the cut and the autoencoder can probe %-level branching ratios.
In the Appendix we show sensitivity projections assuming a much larger 10% systematic background uncertainty, with only very modest degradation in reach. This shows that the overwhelming QCD background has been reduced to low enough levels to make the search statistics limited, and speaks to the robustness of our results.

Conclusion
SUEPs represent a highly plausible but extremely challenging experimental signature of confining hidden sectors, which typically results in a high multiplicity of soft SM final states. To date, there are no targeted LHC searches for SUEP, and most existing searches have limited or no sensitivity. Furthermore, since SUEP is produced by hidden sectors featuring fairly strongly-coupled, approximately conformal dynamics with a wide variety of possible dark hadronization scenarios, modelling the detailed production of SUEP is fraught with uncertainties. Existing proposals for SUEP searches [8,37] target conspicuous, qualitative features of the final state, such as displaced vertices or leptons. Their existence is a prediction for some decay portals of the dark hadrons, and targeting them with searches is fairly robust with respect to modelling uncertainties. An important but previously unaddressed question is how well we can look for SUEP using only its basic kinematic features, without any conspicuous SM final states. This is not only an experimental challenge, the design of such a search also has to be very mindful of uncertainties in the signal modelling. We investigated this worst-case scenario by studying SUEPs from dark hadrons produced in exotic Higgs decays, which decay promptly and purely hadronically. Our choice of model is motivated by the Higgs portal being one of the most plausible production modes for new physics. The modest energy scale of the Higgs decay also eliminates high-energy or high-mass observables as discriminants. Finally, Higgs production in association with a W/Z circumvents the problem of triggering on the SUEP final state directly [37].
Our first results are observables which capture the essential SUEP features from the track momenta. We focused on the charged particle multiplicity N charged , the event ringisotropy [38] I, and the interparticle distance ∆R ij . All three are robust with respect to modelling uncertainties of the SUEP final state, since they capture the essential model features: high multiplicity of final states, and dark hadron production that is more isotropic and democratic in momentum than the QCD background. The charged track interparticle distance matrix ∆R ij is particularly suitable for ML applications.
Based on these observables, we devised three strategies for h → SUEP searches. The first is based on a simple cut on the interparticle distance matrix. The second assumes that our naive signal simulation tools can be trusted, and uses supervised ML techniques. The third is an unsupervised ML approach using a fully connected autoencoder trained only on simulated QCD events as an anomaly detector. Both ML approaches used the slightly modified interparticle distance matrix of Eq. (4.2) as the event representation. The cutand-count approach serves as a simple baseline, over which the unsupervised ML approach represents a significant improvement, even without detailed knowledge of the signal. This is to be compared to the higher sensitivity of the supervised machine learning approach, which is unlikely to be robust with respect to modeling uncertainties or choice of training parameters.
All three approaches will probe exotic Higgs decays to prompt, hadronic SUEP at the HL-LHC for branching fractions at the percent level, with both the autoencoder and supervised approaches probing rates as low as 1%. We assumed a systematic background uncertainty of one percent, but increasing this to ten percent only modestly decreases sensitivity, signaling that our analyses have sufficient differentiation power to reduce the enormous QCD background to a statistics-dominated level. These estimates are conservative, since our simulated background samples were still too small to fully cover the exceedingly large background rejection. A realistic search combining the autoencoder with data-driven background estimates will achieve significantly higher sensitivities.
Our results show that even without a detailed theoretical description of the SUEP showering process, an analysis using an unsupervised neural network can be highly sensitive to exotic Higgs decays to SUEP. Our observables and techniques can equally well be used in SUEP searches using leptons, photons or displaced vertices, to significantly enhance their sensitivity based on the inherently SUEP-y kinematics. Generalizations of our methods, for instance including searches for explicit dark hadron resonances, should allow for sensitivity to SUEPs with higher dark hadron masses or dark Hagedorn temperatures than those targeted by our analysis. Our new, unsupervised search strategy can be applied to a wide range of LHC scenarios to discover new physics, even if the true BSM model should differ from our exact theory expectations. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium [101,102]. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund -Research Excellence; and the University of Toronto. This research was enabled in part by support provided by Compute Canada (www.computecanada.ca). Figure 8 shows the two-dimensional correlations between each pair of the observables N charged , I, and ∆R. While the average values of N charged and ∆R vary in a similar fashion with m D and T D , their distributions for a fixed choice of signal parameters are not highly correlated. Figure 9 shows the AUC for each of the twelve supervised graph networks trained on a single sample of signal events, with parameters indicated by the red dot in each plot. This vividly demonstrates how the the parameter choice of the signal dataset on which the supervised network is trained has a large impact on the regime of parameter space where it can achieve high AUC, providing a strong argument for using the cocktail approach.

A Additional Results
The sensitivity projection for these supervised networks with one percent background systematic is shown in Fig. 10. The peak sensitivity for each network is also at the 1%-level, similar to the other approaches, but for potentially a different range of SUEP parameters, depending greatly on the training parameters of each network. Figures 11 and 12 show sensitivity projections for 10% systematic background uncertainty.