Disorder and Mimesis at Hadron Colliders

We discuss how systems with a large number of degrees of freedom and disorder in their mass matrix can play a role in particle physics. We derive results on their mass spectra using, where applicable, QFT techniques. We study concrete realizations of these scenarios in the context of the LHC and HL-LHC, showing that collider events with a large number of soft b-quark jets can be common. Such final states can hide these models from current searches at the LHC. This motivates the ongoing effort aimed at lowering trigger thresholds and expanding data scouting.


Introduction
The current progress of the experimental effort at the Large Hadron Collider (LHC) has largely exceeded expectations. An unprecedented amount of high-quality data has been collected and true milestones have been reached for the field, such as the discovery of the Higgs boson [1,2]. After ten years of operation, however, and hundreds of measurements that constrain all the most plausible scenarios for physics beyond the Standard Model (SM), we are left to wonder if there are any new particles hiding at the weak scale. This is the right time to ask this question because we have already a rather extensive picture of physics around a TeV from the studies performed on tens of fb −1 of data at √ s = 13 TeV [3, 4]. The upcoming years of LHC operation will be characterized by a very different pace, determined by a slow increase in sensitivity driven by the collected integrated luminosity.
This question has inspired a large part of the theoretical effort for the past several years, mainly in the direction of finding new solutions to the hierarchy problem hidden from traditional LHC searches. By now most of these new scenarios are significantly constrained. Here we would like to take a completely different perspective and abandon all prior theoretical expectations. If we do so, two possible answers strike us as the simplest and potentially most plausible: either we have found nothing beyond the Higgs boson because (1) there is nothing to be found at the weak scale 1 or (2) there are too many new particles.
In this paper we show how large N sectors are naturally hard to detect at hadron colliders. The reasons are simple and independent of a specific model. The first one is that they require a small coupling to be consistent. This together with the finite kinematical range accessible to colliders can give a small total production rate. The second reason is that high-multiplicity final states containing only low p T particles can easily dominate their total production rate, presenting a challenge for current triggers. We discuss both of these aspects in more detail in Sect. 2.
Emphasizing this general point about large N sectors is useful both from the experimental and theoretical perspective. On the experimental side, anything that can be lost due to current triggers deserves serious consideration. Missing new physics because of our own choices in the selection of events would be a highly tragic mistake. An extensive effort in this direction is already in progress at the LHC, in the form of data parking and data scouting [5][6][7][8] and upgrades for the high-luminosity run, such as implementing particle flow at Level 1 [9,10] in CMS and the Feature Extractors [11,12] and Fast Track Trigger (FTK) [13] in ATLAS.
On the theory side, even if the mimetic properties of large N sectors that we discuss in Sect. 2 were not inspired by any open theoretical question, the models that realize them are well-motivated. They can arise as remnants of string theory compactifications [14,15] and/or as a low energy sector of the landscape [16]. They can also be part of a hidden sector containing dark matter or modifying the electroweak phase transition, giving rise to a phase of symmetry non-restoration [17][18][19]. They are related to the broader framework of hidden valleys [20][21][22][23] and realize a phenomenology that is in-between that of traditional hidden valleys (a small number of possibly displaced SM particles in the final state) and that of conformal hidden valleys [24,25] (many particles emitted isotropically).
Furthermore, from a bottom-up viewpoint, the large N sectors that we discuss offer the perfect opportunity to study the possible role of disorder in model building. In our construction disorder is nothing more than a useful phenomenological tool. It allows us to capture possible O(1) variations of the low energy parameters of the theory. However, as it will emerge in the following, disordered systems posses interesting structural properties and in the future their significance in beyond the SM (BSM) physics might be much greater than this. Therefore we use this opportunity to discuss a number of results on random mass matrices in a heuristic way, useful for model building. In the appendices we expand our derivations, making them more rigorous and using, where applicable, path integral techniques familiar from QFT.
The paper is organized as follows: in Sect. 2 we discuss the general kinematical mechanism that hides large N sectors from detection. In Sect. 3 we present a concrete . . . N states are kinematically accessible to a collider gives a production cross section that is suppressed by powers of N LHC /N . realization of the ideas discussed in Sect. 2, in the form of a disordered, large N model of scalars. In Sect. 4 we discuss general results on disordered mass matrices of scalars, in Sect. 5 we study the collider phenomenology of these models, and in Sect. 6 we discuss large N models containing fermions. We conclude in Sect. 7 suggesting next steps to further develop this framework.

Large N Mimesis
The reasons that make a new sector with a large number of new particles hard to detect at hadron colliders are very simple and completely general. We present them briefly in this section before discussing a specific model. Here we imagine that the number N of particles is always large, but the kinematic arguments still apply to a moderate number of particles, even as few as 5.
Introducing N 1 particles in a finite mass range has three main consequences relevant for colliders: 1. The spectrum is compressed and final states with only soft particles are common.
2. The theory requires a small coupling to be consistent.
3. Long decay chains can arise naturally and dominate the total production rate.
The first two items tend to make new sectors with a relatively large number of new particles hard to detect. The last one can be used as a handle to disentangle these new sectors from the background and even reconstruct their structure. However, long decay chains and soft final states go hand-in-hand so having more particles in the final state is not necessarily advantageous. Let us now discuss each of these three aspects in more detail.
Compression is a trivial consequence of having a large number of new states in a finite mass range. However the amount of compression depends only to the density of particles per unit mass, so it can persist also with a small number of new particles. mass small mass long decay chains Figure 2: A theory with a moderately large number of new particles in a finite mass range, having a small coupling to the SM, is characterized by a few light states that decay directly to the SM and events with high multiplicities produced by the long decay chains of the heavier states.
The small coupling arises if we require these theories to be perturbatively consistent. Diagrams that are typically higher-order in the couplings, like loop diagrams, involve sums over the new particles and when their number N is large the couplings must compensate by scaling with the appropriate power of 1/N . In this regime it is useful to reorganize the perturbative expansion as an expansion in powers of 1/N [26].
The small coupling on its own is not enough to guarantee a small production rate. Typically, cross sections scale with the 't Hooft coupling which means that any individual final state is N -suppressed, but the sum over all states can result in an O(1) rate. One exception is when only a subset of the new particles, N LHC , are kinematicallyaccessible. In this case total rates can be suppressed by powers of N LHC /N , as illustrated in Fig. 1.
This leaves us with the last and most interesting aspect of these new sectors: long decay chains. If the new sector is somewhat secluded, i.e. the couplings to the SM are smaller than the couplings between the new states, it is likely for states towards the top of the spectrum to cascade within their own sector before decaying to the SM. In a dense spectrum there is no phase space suppression for not decaying directly to the bottom of the spectrum.
Furthermore, in dense spectra the production rate of states higher in the spectrum can be comparable to that of the lowest-lying states. Since the only guaranteed low multiplicity final states come from production the lightest particle, in such spectra the probability is O(1/N LHC ) for single production and O(1/N 2 LHC ) for pair production so that higher multiplicity final states can dominate the total production cross section.
Observing states near the bottom or top of the spectrum each have their own challenges. If we produce a state towards the bottom of the spectrum, most of the time it decays directly to the SM, generating a low-multiplicity final state already targeted by current searches. However these particles might just be too light to be detected either because of trigger thresholds or because of backgrounds, especially if they decay to jets. Furthermore, the total rate for these low multiplicity final states might just be too small to be detectable. On the contrary a particle near the top of the spectrum cascades within its sector before decaying back to a large number of SM particles. The total visible and/or invisible energy can easily be too small to fire an H T or missing energy trigger if all the new particles are relatively light (i.e. few hundreds of GeV rather than a few TeV) and/or the spectrum is compressed. This is summarized schematically in Fig. 2.
In the next section we introduce a model that makes this simple discussion more concrete and in Sect. 5 we study the statements made in the previous paragraphs quantitatively. We consider sectors with masses around a few hundreds of GeV with N between 5 and 50, but obviously the statements made here are much more general.

A Concrete Model
In this section we introduce an explicit model. Other models can be constructed, for example the one in Sect. 6, but this one is the simplest and illustrates all relevant features. Consider N real scalars with the Lagrangian where sums over repeated indices are implied and run from 1 to N . The parametrization should be interpreted as an expansion around a local minimum, not necessarily valid for arbitrarily large field excursions. In addition to Eq. (1), we connect the new scalars to the SM through the most relevant interactions that are possible, the Higgs portal couplings, where again sums over repeated indices are left implicit. In Sect. 5, where we discuss collider phenomenology, we consider models with N ranging from 5 to 50. We study separately three distinct phenomenological possibilities: 2. All interactions are present and single production dominates.
3. All interactions are present and pair production dominates.
The first case has a potential that respects a Z 2 symmetry under which It can be considered our "nightmare scenario" being maximally difficult to detect at colliders. The second one is a more faithful description of a friendly landscape [16], where the new scalars can get vacuum expectation values of the order of their masses. The third possibility covers a different limit of this model that has distinct phenomenological features. From the point of view of the Lagrangian it is similar to the first case, but with the addition of a small Z 2 -breaking. The breaking is sufficiently small that pair production is still the dominant production process at the LHC. The model defined by Eqs. (1) and (2), maps onto models used for baryogenesis (with a particular choice of interactions) [17][18][19] and can be a QFT model of the landscape [16,[27][28][29][30]. The new scalars can also be moduli from extra dimensions compactified at some large scale M * . If the compactification scale is much larger than their mass, they would have to be tuned moduli to be visible at colliders, i.e. they need couplings not suppressed by some power of m φ /M * . Given the ubiquitous presence of the weak scale in nature their presence around LHC energies might not be a coincidence. The dark matter energy density today is ρ DM ∼ (v 2 /M Pl ) 4 , the cosmological constant and neutrino masses are also related to the same combination of v and M Pl : Λ 1/4 CC ∼ m ν ∼ v 2 /M Pl , not to mention the role of the weak scale in the SM itself. From a more pragmatic perspective, this might be just one of many sectors spread over many orders of magnitude in mass that arise from the compactification of extra-dimensions and supersymmetry breaking.

Large N
To have a well-behaved perturbative theory the couplings must scale as some inverse power of N if N 1. For the Z 2 -symmetric case with a ijk = a SM i = 0, at large N we need λ SM 4π/N and λ 16π 2 /N 2 . This can be verified by inspecting diagrams. For the kind of arguments needed to derive this scaling in general see Ref. [31]. When trilinear couplings are included we additionally require a SM 4πv/ √ N and a 4πm φ /N 3/2 , where m φ is the typical mass scale of the scalars and v ≈ 174 GeV the vacuum expectation value of the SM Higgs.

Randomness
We introduce phenomenological randomness into the model through the mass matrix. As discussed in the next section we choose fully-populated matrices where each entry is drawn randomly from the same Gaussian distribution. For N real scalars the mass matrix is then symmetrized. In Eq. (1) we wrote down the model after mass diagonalization. When considering the phenomenology of the model in Sect. 5 we take advantage of the fact that in the large N limit there is a known distribution for the mass eigenvalues. For simplicity, we start from the diagonal basis and draw the masses randomly from the known spectral density discussed in the next section.
We take the trilinear and the quartic couplings to be roughly of the same order. This reflects our choice of a fully-populated mass matrix with entries drawn from the same distribution, as in this case all mass eigenstates overlap at O(1/ √ N ) with all flavor eigenstates.

Existing Constraints
The most stringent bound on these models is indirect and arises from LHC measurements of Higgs couplings. The mixing between the Higgs and the new scalars reduces all Higgs couplings, modifying the global signal strength measured at the LHC. We can estimate the impact of this bound by taking all scalar masses to be approximately m φ . Then we find the modification of Higgs couplings to be governed by the parameter and are constrained at permille level making them subdominant [36]. Above we have again imagined all scalars to have the same mass m φ .
The new states can also be singly produced at colliders through the mixing with the Higgs. The most relevant direct searches are those for heavy Higgses, but existing ones, once the bound from Higgs coupling measurements is taken into account, have no residual exclusion power [37][38][39][40][41]. However, the full 13 TeV dataset has not yet been analyzed and resonant searches for pairs of electroweak gauge bosons will become relevant in the future.
In the Z 2 -symmetric case the new scalars do not mix with the Higgs boson and the constraints that we have just discussed do not apply. Furthermore, the new particles can only be pair produced at colliders with much smaller rates. Current searches for pairs of new particles have not even begun to probe λ SM ∼ O(1) [42][43][44][45].

Disordered Model Building
In this section we discuss the spectrum and eigenvectors of disordered mass matrices. We focus on N × N matrices with entries randomly drawn from Gaussian distributions. Our starting point is the Gaussian Orthogonal Ensemble (GOE), defined by where we use M ii to denote the diagonal entries of the mass matrix M and M ij for the off-diagonal entries. As shown in Eq. (6), the entries are drawn from the unit Gaussian N (0, 1) with mean µ = 0 and standard deviation σ = 1. 2 The GOE contains real symmetric matrices and thus models of real scalars. For complex scalars the analog ensemble is the Gaussian Unitary Ensemble (GUE) that contains complex Hermitian matrices. Most of our results apply to both cases, with differences that we specify in the following. The natural generalization of the GOE, with broader physical applications, that we consider in this section is where the diagonal entries are now drawn from N (µ d , σ d ) and the off-diagonal entries from N (µ o , σ o ). In the next subsection we often take µ d = µ o = µ.

Summary of Main Results
The main result of this section is that in the large N limit the spectral density of mass matrices with Gaussian-random entries follows a universal distribution called the Wigner semicircle distribution [46][47][48], which once appropriately normalized reads For finite N and including the Gaussian parameters of Eq. (7) (but fixing We have omitted from this equation a single large eigenvalue located at N µ. The value of β depends on the ensemble: β = 1 for the GOE and β = 2 for the GUE. To obtain Eq. (8) from Eq. (9) we have rescaled σ o in the large N limit. In Fig. 3 we show a comparison of the Wigner semicircle distribution with the spectrum of two GOE matrices. The eigenvectors of these matrices do not significantly deviate from what one might naively expect. Since the mass matrix is fully populated with elements of comparable magnitude, we can imagine the eigenvectors to be an approximately democratic mix of flavor eigenstates with weight ∼ 1/ √ N . This can be confirmed mathematically [49]. For example the squares of the components of the GOE eigenvectors follow a beta distribution with mean 1/N and variance of O(1/N 2 ) [49]. There is not much more to say regarding eigenvectors that we found useful for phenomenology. On the contrary Wigner's semicircle law has several interesting consequences: 2 The standard deviation of the off-diagonal entries is 1/ √ 2 rather than 1 because we imagine that the matrix is symmetrized via (M + M T )/2. • All eigenvalues except the largest one are roughly contained in the interval |m| √ N σ o . In particular if σ o µ they can be much lighter than the scale µ of the matrix components.
• The typical eigenvalue spacing is of order σ o / √ N .
• The lightest eigenvalue has a mass of order The last two are trivial consequences of having N eigenvalues in a 2 √ N σ o interval centered around zero. Instead the first observation is more interesting and seems extremely useful for model building. We can imagine that µ ∼ M Pl , while σ o is much smaller, maybe around the weak scale. This naturally generates a hierarchy. However we have only rewritten N −1 Goldstone bosons in an unusual basis. Take N real scalars, contained in a vector Φ, with an O(N )-symmetric potential If we expand V (Φ) around the true minimum and take the vev in the direction Φ ∝ This explains why the typical scale of the eigenvalues is σ o which is explicitly breaking the symmetry, instead of µ that preserves it. The reason why σ d does not play a role at large N is simply that we have N diagonal elements and ∼ N 2 off-diagonal ones.
We will discuss this more precisely in the next subsection. The last useful result on the spectrum is that taking a different mean for the diagonal elements µ d , shifts the spectral density from zero by the amount (µ d − µ): If viewed in terms of Goldstone bosons this is a simple consequence of the explicit breaking of the O(N ) symmetry. Before turning to the derivation of these results it is worth to comment on a difference between our models and disordered condensed matter systems. It is natural to ask if we can pick any variance for our mass matrix or if we need something scaling as an inverse power of N to have a consistent theory. In real disordered systems locality and the central limit theorem imply that the free-energy is self-averaging, i.e.
where · is a disorder average. The argument is quite intuitive: only particles that are nearby interact with each other appreciably and we can divide F into a sum over many cells that are not strongly interacting with each other. In the large N limit we can ignore surface effects and recover the result above. If this result applied to us, it would require for example σ 2 o ∼ 1/N for our scalars, just from computing the free energy from the partition function Z However in our theories there is no notion of locality in the same sense as for condensed matter systems. The scalar with flavor 1 can mix as strongly with that of flavor 2 as with any other, so in this sense all interactions are long range.

Catalan Numbers
Before turning to the derivation of the above results we find interesting to discuss one structural property of our large N sectors with Gaussian mass matrices. The set of numbers known as the Catalan numbers, impacts both the shape of the eigenvalue distribution and the length of our decay chains. The moments of Wigner's semicircle distribution are the Catalan numbers A derivation of this result can be found in App. B while an explanation that uses a remarkable connection with planar diagrams is presented in App. C. To see the role of the Catalan numbers for our decay chains, consider a scalar sector with particles that can either decay to two other scalars or to two SM particles. This is a good approximation of our general model in Sect. 3 since we expect the trilinear couplings to dominate the branching ratios.
After a first decay, each daughter scalar would then likewise either decay to additional scalars or to two SM particles. Each possible final state can then be represented by a binary tree. The average number of final state particles can be approximated by computing the weighted sum of possible final states. This requires knowing, for a given number of leaves 2n, how many distinct binary trees there are. This sequence of numbers is again given by C n .
The asymptotic behavior of the Catalan numbers also lets us make a rudimentary estimate of the average decay chain length. Let the probability of decaying to the SM be p and the probability of decaying to the new sector be q. The average number of final state particles, ignoring all phase space factors, is then Nmax n=0 (n + 1)C n (pq) n p (16) with N max determined by phase space. For n → ∞ we have that C n ∼ 4 n /n 3/2 which means that the nth term in the average goes like ∼ (4pq) n / √ n. Assuming that (4pq) ∼ O(1) the average will go like ∼ √ N max . N max is set by phase space and is O(100) for the parameters considered in Sect. 5 in good agreement with the numerical results in the same section. This heuristic derivation is valid for two-body decays, but the scaling of ∼ √ N max continues to apply for higher n-body decays.

The Joint Eigenvalue Distribution
Most of the results in the previous section can be derived using path integral techniques and can be estimated by simple dimensional analysis. In this section we go through the heuristic arguments that justify the form of the spectrum and point to the appendices where a more complete QFT-inspired derivation is presented. The first step is to go from the distribution of matrix components in Eq. (7) to the joint distribution of eigenvalues ρ(m 1 , ..., m N ) via the change of basis where ds 2 M is the distance on the space of matrices, more details are given in App. A. The other differentials are The joint distribution of eigenvalues is then ρ(m 1 , ..., m N ) = ρ M (M D )|J(M D )| dU and we would like to have a closed-form expression for it. If we think about ρ M as an action for a matrix model, the eigenvectors U represent the gauge freedom associated to the choice of basis and J is a gauge-invariant Fadeev-Popov determinant. Then it is not surprising that J does not depend on U . The proofs of these statements can be found in App. A.
If J depends only on the eigenvalues we can make an ansatz for its form based on the following arguments. When two eigenvalues coincide the transformation becomes singular and J must be zero, so we expect To determine β we can use dimensional analysis 3 In App. A we compute J explicitly. Now we are left with evaluating ρ M (M D ). If we take µ d = µ o = 0 in Eq. (7) it is easy to conclude that in the large N limit 4 Combining this with the result for the Jacobian we finally obtain the joint distribution of eigenvalues This result becomes exact in the limit σ o = σ d / √ 2 [50]. From Eq. (22) is possible to derive the Wigner semicircle distribution, Eq. (8), by solving a path integral [50,51]. We can either use Feynman diagrams, as done in App. C, or use a saddle point approximation [46,48]. Even without going through all the derivation we can understand most of the results in the previous section just by looking at Eq.
and taking the derivative dρ/dm * to find the maximum. As pointed out before the other relevant statements on the spectrum descend trivially from having N eigenvalues in an interval that is roughly 2 √ N σ o wide. We derived these results by assuming µ d = µ o = 0, but it is easy to show that they are valid also for µ d = µ o ≡ µ = 0. In this case we can split M into a zero-mean component M and a matrix A with all equal entries, If we apply an orthogonal (or unitary depending on which ensemble we are considering) transformation U A that diagonalizes A to the right hand side of the previous equation, we are left with a matrix of zeros plus one large eigenvalue N µ: By diagonalizing A, M is rotated as before, while the term proportional to the identity, (µ d − µ o )1 N ×N , is unaffected. Therefore we still have Wigner's semicircle law, but shifted from zero by an amount (µ d − µ o ). In addition we have the large eigenvalue which is also shifted by a small amount m max = N (µ o − 1) + µ d . This concludes our heuristic derivation of results on the spectrum of GOE and GUE matrices. We now turn to another interesting phenomenological possibility: having mass matrices that mix only p nearest neighbors.

Band Matrices
It is interesting to consider what happens if we draw the elements of our matrices from the same probability distributions considered in the previous sections, but instead of populating the full matrix we allow only for nearest neighbor interactions among p < N 11 m 2 12 · · · m 1p · · · 0 0 m 2 12 m 2 22 · · · · · · · · · 0 0 . . · · · · · · · · · · · · · · · · · · m 2 (N −1)N 0 0 · · · · · · · · · m 2 This scenario can arise from localization in an extra dimension for example. The eigenvalues of these matrices are spread over a smaller range than those in the GOE or GUE. ρ M (M ) is the same as before (see Eq. (7)), leaving intact the Gaussian measure in Eq. (22), but the Jacobian of the trasformation, even just on dimensional grounds, cannot be the same. If we repeat the arguments in the previous subsection we expect β to be wide. This is indeed confirmed numerically, as shown in Fig. 4. Also the eigenvectors of these matrices, not surprisingly, are more localized than their GOE or GUE counterparts. The limiting case is a diagonal matrix that exhibits perfect localization. In the opposite limit we might have eigenvectors as spread as in the previous cases. If the off-diagonal elements are smaller than the diagonal ones we can have a stronger form of localization, known as Anderson localization [52], with mass eigenstates spanning approximately p flavor eigenstates. This fact has already found several applications in the context of BSM physics [53,54] and cosmology [55][56][57].
We do not explore in detail band matrices in this paper, but from the discussion in this section we expect a collider phenomenology similar to the one that we discuss in Sect. 5 with potentially longer decay chains compared to fully populated matrices.

Collider Phenomenology
In this section we explore the collider signals predicted by the scalar model that we presented in Eqs. (1) and (2). In different regions of parameter space, quite different behavior is expected, ranging from simple dijet resonances (singly or pair produced) to long cascades ending in many SM particles and possibly missing energy. In all cases presented here we extract particle masses from a Wigner semicircle distribution that has support between 100 and 600 GeV.
We study three regions of parameter space that provide a good representation of the possible final states of the model. The first case is the Z 2 -symmetric theory where the only coupling to the SM is a quartic Higgs portal coupling, λ SM . The scalars are pair produced through this coupling (Fig. 5a) and can decay either through the same Higgs portal coupling (Fig. 6b) or through the hidden sector quartic λ (Fig. 6a). For simplicity we take λ SM = 1/N and λ = 1/N 2 for all the scalars 6 .
In the second case we allow non-zero trilinear terms to be present. This changes both the production and the dominant decay channels of the scalars. Production occurs both through pair production (Fig. 5a) and single production (Fig. 5b). The scalars can decay either as before ( Fig. 6a and Fig. 6b), or through the two-body scalar channel (Fig. 6c), or to a pair of SM particles via mixing with the Higgs (Fig. 6d). For the quartic couplings we use λ SM = 0.1/N and λ = 1/N 2 . For the trilinears we take a = m min /N 3/2 , where m min is the lightest scalar mass, and a SM = a, for all the scalars.
In the last scenario we still allow all the couplings to be present, but we take very small trilinears: a = 10 −5 × m min and a SM = 10 −5 × m h , leaving the quartics as in the exactly Z 2 -symmetric case. In this case pair production dominates over single production and the lightest scalar decays to a pair of SM particles.
Before discussing the phenomenology of these three scenarios it is useful to take a look at the scalar production rate at the LHC. We show the total event rate from gluon fusion summed over all pairs of scalars in Fig. 7 (left) as a fuction of the coupling λ SM and summed over all singly produced scalars in Fig. 7 (right) as a function of the trilinear a SM . We use a center-of-mass energy of 14 TeV and an integrated luminosity of 3 ab −1 . Cross sections were calculated using Madgraph 5 [58] with gluon fusion implemented via the Higgs Effective Theory module [59]. We have plotted Fig. 7 for a representative choice of the scalar spectrum drawn at random from a Wigner semicircle distribution in the range of 100 GeV to 600 GeV. The red lines are for N = 10 and the purple lines are for N = 50. The figure allows us to conclude that even with our N -suppressed couplings we still have a reasonable number of events to work with at the end of the high-luminosity program of the LHC.
To analyze the phenomenological features of the model we start by showing the total number of particles N tot per event. In Fig. 8 we show N tot for N = 5, N = 10, and N = 50 scalars in the Z 2 -symmetric case. In Fig. 9 we show the same for the second case, where single production dominates, while in Fig. 10 we show N tot for the last scenario that is nearly Z 2 -symmetric, but has small trilinears that allow the lightest scalar to decay. The different colors correspond to a different randomly-generated mass spectra, so the four histograms in each figure represent the variation that we expect from randomness in the mass matrix. Especially at small N , the variation between different spectra is primarily a function of the gap between the lightest state φ N and the states closest in mass to it. When there is a large gap between φ N and the next state, events with N tot = 2 are favored because φ N φ N production dominates the overall rate. Spectra for which the N tot = 2 bin is small have states that are near in mass to φ N so that their production rate is comparable to that of φ N . Furthremore, smaller mass gaps near the bottom of the spectrum mean that heavy scalars have a very similar probability to decay to any of the light scalars favoring longer decay chains. If φ N is much lighter, phase space will instead favor a direct decay to the bottom of the spectrum.
This discussion leads us to identify the main qualitative feature common to all three regions of parameter space: increasing the density of states in a fixed mass interval, long decay chains become more common, giving rise to higher multiplicity final states. This emerges clearly from Fig. 11 where we compare one of the spectra of Fig. 8 for N = 50 with N = 5 scalars in the same scenario, but distributed over a much smaller mass range: from 200 to 260 GeV. The two spectra have the same average mass splitting between neighboring states and very similar final state multiplicities.
From Figs. 8, 9, and 10 another general aspect of the phenomenology of the model emerges clearly: when single production dominates traditional resonant searches for    tion dominates the main signature consists in four or more particles in the final state. In the Z 2 symmetric case these particles are mostly soft b-quark jets. This is shown in Fig. 12 and Fig. 14. The first figure counts the number of b-quarks in the event, while the second one shows the fraction of events with average energy per particle above a certain threshold for N = 10 and N = 50 scalars. From Fig. 14 it is clear that our final states are extremely challenging for traditional low-multiplicity triggers and even the total energy in the event, shown in Fig. 15 is not a good handle. A quantitative discussion of trigger thresholds goes beyond the scope of this work, but it is clear from our results that low-threshold, high-multiplicity triggers are well motivated in this scenario. Note also that in this scenario the lightest scalar φ N is stable and we have always a small amount of missing energy in the event. It is too small for triggering purposes, but it can be used as a handle to identify this model.  Figure 14: Fraction of events with average energy per particle above E cut . The plots were made in the Z 2 symmetric scenario where a SM = a = 0 in Eqs. (1) and (2). Left: 10 scalars, Right: 50 scalars.
Let us now turn to the last scenario, where pair production dominates but φ N can decay to a pair of SM particles. Interestingly for N 10 most of the events contain at least four W 's, as shown in Fig. 13. At larger N b-quarks are still the dominant species of SM particles. In this scenario φ N decays via its mixing with the Higgs. So the large W multiplicity is due to the mass of the lightest state in the spectrum that for small N is usually larger than the Higgs mass (due to the shape of Wigner's semicircle). The W W decay width turns on rapidly for m W m φ 2m W . On the contrary at larger N when we start to populate the entire Wigner distribution, we always find a scalar with mass comparable or smaller than m h for which decays to SM quarks dominate. From the point of view of triggering the total energy in the event still does not offer a very good handle as shown in Fig. 17. However at small N the average energy per particle (shown in Fig. 16) is more than enough for leptonic and some multijet triggers.
To summarize, the phenomenology of the model is very rich, ranging from scenarios where traditional resonant searches capture the bulk of the events to cases where long decay chains with multiple b's or W 's are the most common signatures. In general the total energy and missing energy in the event can not be used for triggering and high-multiplicity triggers are motivated. In this section we took the mass range of the scalars between 100 and 600 GeV. It would be interesting to explore different mass ranges. Going to larger masses would boost the total energy in the event potentially changing our qualitative conclusions on triggering. However this can be done only at the price of considerably reducing the production rate and would not make these new sectors less elusive. Going to lower masses would make the new scalars even harder to trigger on or move them to kinematical regions better explored outside of the domain of hadron colliders. However following the agnostic approach outlined in the introduction it would be worth to consider also much lighter sectors and a completely different set of experiments. We leave this exploration to future work.

Fermionic Hidden Sectors
In this paper we have chosen to focus on models with scalars and explore thoroughly their mass matrices. However the same ideas could be realized in models containing a large number of fermions or a mixture of particles with different spins. In this section we discuss some of the differences that one would encounter for fermions. The chiral protection of fermion masses makes them plausible low energy remnants of the compactification of extra dimensions, also if the typical scale of compactification is much larger than their mass. For example low energy fermions are common in some realizations of string theory [14]. One example that would have a collider phenomenology similar to the theories in Sect. 3 is a model of N Weyl fermions with a Yukawa interaction This choice is rather appealing from the point of view of coupling this sector to the SM since it gives us a clean symmetry argument to forbid the lHψ α vertex leaving us with the interaction which gives the long decay chains discussed in Sect. 5. If we also take the mass of φ in a range that makes three body decays within the dark sector (ψ α → ψ β φ * → ψ β ψ γ ψ δ ) relevant we can reproduce exactly the same structure that we had in the scalar models.
However when considering the mass matrices of N fermions one should be aware of some differences compared to the scalar case. We have to diagonalize M † M to obtain the pole masses in the free theory, while the parameters in the Lagrangian form a linear mass matrix M . If we treat the the entries of M as the fundamental parameters of the theory, drawn from a random distribution, the asymptotic form of the spectral density is different from what discussed for GOE and GUE.
Given a square N ×N matrix M , with i.i.d. entries drawn from a Gaussian distribution with zero mean and unitary variance the asymptotic distribution for the eigenvalues of (M † M )/N is a particular case of the Marčenko-Pastur density [60] which denotes an accumulation of eigenvalues around zero. A diagrammatic derivation of this result can be found in [61]. Note that we have used m for the eigenvalues of M † M that have the dimension of a squared mass. The Marčenko-Pastur density can be easily generalized to the case of non-unitary variance. Assuming that all the entries of M have zero mean and variance σ, the asymptotic distribution for the spectral density of (M † M )/N is Also in this case turning-on a non-zero mean µ equal for all the entries does not affect the spectral density except for the appearance of one large eigenvalue of M that is O(N µ). We leave a more detailed discussion to future work.

Outlook
In this paper, motivated by the current null results at colliders, we have ignored some of the unspoken rules of BSM model building.
We have asked what is the simplest scenario in which new particles are present at the weak scale, but still invisible at colliders, without attempting to answer the open questions that have driven the field. Even if our starting point was orthogonal to most traditional phenomenological studies, our "kinematics of invisibility" is realized in a number of scenarios that are well motivated theoretically.
More concretely, we have discussed the simple kinematical consequences of having a large number of new particles in a finite mass range. This situation can arise in many BSM scenarios and has characteristic phenomenological signatures that we have explored in Sect. 5. The main messages are that final states with multiple soft particles can be common and that the total production cross section can easily be small enough to motivate the high luminosity upgrade of the LHC. Often there is not enough total energy to pass global triggers such as those targeting high H T or missing E T signatures. This provides further motivation for the ongoing work aimed at lowering trigger thresholds.
The new sectors that we have studied are naturally described in terms of disordered mass matrices and couplings. This is a convenient phenomenological tool to parametrize our ignorance. It allows us to keep parameters that are supposed to be of the same order close to each other without giving up potentially interesting O(1) variations. The main phenomenological interest lies in accidental compressions near the bottom of the spectrum. We find that this situation is not uncommon and can considerably increase the number of soft final state particles in the event. Aside from this point, disordered mass matrices have an interesting structure that we have discussed in Sect. 4 and expanded upon in the appendices. We have made an effort to rederive all relevant results in a language as close as possible to QFT. We hope that further explorations of disorder in model building will have interesting implications for the long standing questions in the field. Even if just at the level of intriguing coincidences we already find much more structure than we naively expected.
There are a number of new directions that this work suggests. Most of them are simple such as expanding the analysis of fermion models and of the combinatorial properties of long decay chains. However the one that we find most intriguing is the general exploration of large N sectors and disorder in phenomenology, especially their aspects that we have not touched in this paper as the possibility of having a large number of metastable vacua and glassy phases.
In conclusion we have presented a simple phenomenon that motivates new explorations of hadron collider data, found connections with motivated BSM scenarios, and introduced some of the tools of large N disordered models in a particle physics context.

A The Vandermonde Determinant
As discussed in Sect. 4 if we take a N ×N Gaussian matrix M , whose matrix components M ij are drawn from the joint probability distribution of its eigenvalues m i has the form where β = 1 for a symmetric matrix and β = 2 for a Hermitian matrix. Z N,β is a normalization factor and the exponential term in Eq. (33) was derived in Sect. 4. The last factor |J| = i>j |m i − m j | β is the Jacobian of the transformation taking us from matrix components to eigenvectors and eigenvalues, it is also known under the name of Vandermonde determinant. It is straightforward to derive it by performing the change of basis To compute J we need to recall the definition of the metric on the space of symmetric matrices where we have used d(OO T ) = 0 → dO T = −O T dOO T , which also shows that O T dO is an antisymmetric matrix. This gives us a metric tensor and we can use the square root of its determinant to obtain J in Eq. (34). Note that when we write ρ M (M )DM we are integrating only over the N (N + 1)/2 independent variables in M . Given the form of ρ M , integrating over the other matrix components would just give an overall constant absorbed by the normalization. So we do not need all the components of the metric tensor defined by Eq. (34) to compute J. With this caveat, combining Eq. (34) and Eq. (35) we obtain It is not hard to generalize these steps to the case of Hermitian matrices.
We find instructive to derive J also as a Fadeev-Popov determinant. The correlation functions of gauge invariant operators O i (A) in a Yang-Mills theory have the same structure as the expectation values of quantities that depend only on the eigenvalues of M where δ (N 2 ) is an N (N + 1)/2-dimensional Dirac delta for symmetric matrices and a

B Moments of Wigner's Semicircle
Recall the Wigner semicircle distribution: The odd moments of this distribution vanish by symmetry. Here we show that the even moments, µ 2n , for integer n, are the Catalan numbers, C n , where C n = 1 n + 1 2n n .

C The Catalan Numbers and Planar Diagrams
Here we present a partial derivation of the Wigner semicircle distribution using Feynman diagram techniques. This derivation highlights the connection between Catalan numbers and planar diagrams.
Recall that the Catalan numbers are the even moments of the Wigner semicircle distribution Consider a Hermitian matrix M . 7 It has an associated Green's function where m i are the N eigenvalues of M . We are interested in the average over many realizations of M The averaging merges the poles into a cut spanning the support of ρ(m). From the Green's function one can find the spectral density ρ(m) via the identity We compute G(z) using the Feynman diagram expansion developed in [62] and discussed pedagogically in [51]. It is convenient to first introduce G i j (z) as We expand G i j (z) to find Odd powers of M vanish in the average. The numerator of each term, explicitly, is where Z is a normalization factor. Eq. (55) resembles a path integral for the matrix M . Computing the full propagator requires evaluating these integrals which we can do using Feynman diagrams. Using Feynman diagrams we can see that in the large N limit planar diagrams dominate, similar to large N QCD [26]. The Feynman rules are shown in Fig. 18. The n = 1 term is shown in Fig. 19. Given that there is no integral over space or time, each diagram contributes a pure number to (M 2n ) i j . Closed loops correspond i j Figure 19: The n = 1 term from Eq. (53) in the computation of G i j (z).
to factors of N from tracing so that non-planar diagrams are suppressed by powers of 1/N . By inspecting diagrams, such as the one in Fig. 19, we can conclude that there are as many planar diagrams with n vertices as there are non-crossing partitions of a lattice with n sites. The definition of a non-crossing partition is precisely that if one puts the n points of a lattice on a circle and connects each point with the next member of its part by an internal path (in cyclic order), the paths do not cross. The Catalan numbers count, among other things, non-crossing partitions [63]. Therefore With Eq. (54) this implies which was shown explicitly in App. B.