Disorder and Mimesis at Hadron Colliders

TL;DR

This paper explores how disordered systems with many degrees of freedom can influence particle physics, particularly at colliders, revealing potential hidden signals in complex final states that current searches might miss.

Contribution

It introduces a framework for analyzing disordered mass matrices in particle physics and applies it to collider scenarios, highlighting new signatures at the LHC and HL-LHC.

Findings

Collider events with many soft b-quark jets are common in these models.

Such events can evade current LHC searches due to their complexity.

Lowering trigger thresholds can improve detection of these hidden signals.

Abstract

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.

Equations

$${\mathcal{L}}_{\varphi }=\frac{1}{2}{\partial }_{\mu }{\varphi }_i{\partial }^{\mu }{\varphi }_i-\frac{m_i^2}{2}{\varphi }_i^2-a_{ijk}{\varphi }_i{\varphi }_j{\varphi }_k-{\lambda }_{ijkl}{\varphi }_i{\varphi }_j{\varphi }_k{\varphi }_l,$$

$${\mathcal{L}}_{\varphi H}=-a_i^{\mathrm{SM}}{\varphi }_i|H{|}^2-{\lambda }_{ij}^{\mathrm{SM}}{\varphi }_i{\varphi }_j|H{|}^2,$$

$${\varphi }_i\to -{\varphi }_i,\forall i.$$

$$\sin \theta \approx \frac{a^{\mathrm{SM}}v\sqrt{N}}{m_{\varphi }^2-m_h^2}.$$

$$\hat{S}\approx \frac{\alpha }{48\pi s_w^2}{\sin }^2\theta \log \frac{m_{\varphi }^2}{m_h^2},\hat{T}\approx -\frac{3\alpha }{16\pi c_w^2}{\sin }^2\theta \log \frac{m_{\varphi }^2}{m_h^2},$$

$$\rho [M]=\prod _{i=1}^N\frac{e^{-(M_{ii}{)}^2/2}}{\sqrt{2\pi }}\prod _{i<j}\frac{e^{-(M_{ij}{)}^2}}{\sqrt{\pi }},$$

$$\rho [M]=\prod _{i=1}^N\frac{e^{-\frac{(M_{ii}-{\mu }_d{)}^2}{2{\sigma }_d^2}}}{\sqrt{2\pi }{\sigma }_d}\prod _{i<j}\frac{e^{-\frac{(M_{ij}-{\mu }_o{)}^2}{2{\sigma }_o^2}}}{\sqrt{2\pi }{\sigma }_o}.$$

$${\rho }_{\mathrm{SC}}(m)=\frac{1}{2\pi }\sqrt{4-m^2}.$$

$$\rho (m)=\frac{1}{\pi \beta N{\sigma }^2}\sqrt{2\beta N{\sigma }^2-m^2}+\mathcal{O}(1/N).$$

$$V(\Phi )=-\frac{m^2}{2}{\Phi }^T\Phi +\frac{\lambda }{4}({\Phi }^T\Phi {)}^2.$$

$$\rho (m)=\frac{1}{\pi \beta N{\sigma }_o^2}\sqrt{2\beta N{\sigma }_o^2-(m-({\mu }_d-\mu ){)}^2}+\mathcal{O}(1/N).$$

$$⟨F^2⟩-⟨F{⟩}^2=\mathcal{O}(1/N),$$

$$F=-T\log Z.$$

$$C_n=\frac{1}{n+1}\left(\genfrac{}{}{0ex}{}{2n}{n}\right),$$

$$C_n={\int }_{-2}^2{\rho }_{\mathrm{SC}}(m)m^{2n}dm.$$

$$\sum _{n=0}^{N_{\max }}(n+1)C_n(pq{)}^{n}p,$$

$$M=U{M}_D{U}^{†},$$

$$\rho [M]DM=\rho [M_D]|J(M_D)|DmDU.$$

$$d{s}_M^2=\mathrm{Tr}[dMdM],$$

$$Dm=\prod _{i}d{m}_i,DU=\prod _{i>j}(U^{†}dU{)}_{ij}.$$

$$|J|=\prod _{i<j}|m_i-m_j{|}^{\beta }.$$

$$\mathrm{GOE}:$$

$$\mathrm{GUE}:$$

$$\rho [M]\propto e^{-\frac{1}{2{\sigma }^2}\mathrm{Tr}[M^2]+\mathcal{O}(1/N)}=e^{-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^N{m}_i^2+\mathcal{O}(1/N)}.$$

$$\hat{\rho }[m_1,\dots ,m_N)]\equiv \rho (M_D)|J(M_D)|=\frac{1}{Z_{N,\beta }}{e}^{-\frac{1}{2{\sigma }_o^2}{\sum }_{i=1}^N{m}_i^2}\prod _{i<j}|m_i-m_j{|}^{\beta },$$

$$\hat{\rho }[m_1,\dots ,m_N]\sim e^{-\frac{1}{2{\sigma }^2}{m}_{\ast }^2}|m_{\ast }{|}^{\beta (N-1)},$$

$$M=M+A-A=M^{\prime }+A.$$

$$M=M+B-B=M^{\prime }+B=M^{\prime }+A+({\mu }_d-{\mu }_o)1_{N\times N}.$$

$$M_{p}^{2}=\left(\begin{array}[]{cccccc}m_{11}^{2}&m_{12}^{2}&\cdots&m_{1p}^{2}&\cdots&0\\ m_{21}^{2}&m_{22}^{2}&\cdots&m_{2p}^{2}&\cdots&0\\ \vdots&\vdots&\ddots&\cdots&\cdots&\vdots\\ m_{p1}^{2}&m_{p2}^{2}&\vdots&\ddots&\cdots&\vdots\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\cdots&\cdots&m_{NN}^{2}\end{array}\right).$$

$${\mathrm{GOE}}_p:$$

$${\beta }_p=\frac{p(p+1-2N)}{N(1-N)}\approx \frac{2p}{N}.$$

Full text

1. 1 Introduction

2. 2 Large $N$ Mimesis

3. 3 A Concrete Model

4. 4 Disordered Model Building

5. 4.1 Summary of Main Results

6. 4.2 The Joint Eigenvalue Distribution

7. 4.3 Band Matrices

8. 5 Collider Phenomenology

9. 6 Fermionic Hidden Sectors

10. 7 Outlook

11. A The Vandermonde Determinant

12. B Moments of Wigner’s Semicircle

13. C The Catalan Numbers and Planar Diagrams

1 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 $\sqrt{s}=13~{}{\rm 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 scale111Including the possibility that new physics is too weakly-coupled for us to detect it. 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 this 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 are 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.222See [26] for uses in dynamical dark matter [27, 28]. In our construction disorder is nothing more than a useful phenomenological tool. It allows us to capture possible $\mathcal{O}(1)$ variations of the low energy parameters of the theory. However, as it will emerge in the following, disordered systems possess 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 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 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.

2 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\gg 1$ particles in a finite mass range has three main consequences relevant for colliders:

The spectrum is compressed and final states with only soft particles are common. 2. 2.

The theory requires a small coupling to be consistent. 3. 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 on the density of particles per unit mass, so it can persist also with a small number of new particles.

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$ [29].

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 $\mathcal{O}(1)$ rate. One exception is when only a subset of the new particles, $N_{\rm LHC}$, are kinematically-accessible. In this case total rates can be suppressed by powers of $N_{\rm 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 the production of the lightest particle, in such spectra the probability is $\mathcal{O}(1/N_{\rm LHC})$ for single production and $\mathcal{O}(1/N_{\rm LHC}^{2})$ 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. a 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.

3 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

[TABLE]

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,

[TABLE]

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:

$a_{ijk}=a^{\rm SM}_{i}=0$. 2. 2.

All interactions are present and single production dominates. 3. 3.

All interactions are present and pair production dominates.

The first case has a potential that respects a $Z_{2}$ symmetry under which

[TABLE]

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 [30, 16, 31, 32, 33]. 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_{\phi}/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 $\rho_{\rm DM}\sim(v^{2}/M_{\rm Pl})^{4}$, the cosmological constant and neutrino masses are also related to the same combination of $v$ and $M_{\rm Pl}$: $\Lambda_{\rm CC}^{1/4}\sim m_{\nu}\sim v^{2}/M_{\rm 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 when $N\gg 1$ the couplings must scale as an inverse power of $N$. First considering only the quartic interactions, we need that $\lambda^{\rm SM}\lesssim 4\pi/N$ and $\lambda\lesssim 16\pi^{2}/N^{2}$. This can be verified by inspecting diagrams. For the kind of arguments needed to derive this scaling in general see Refs. [34, 35].

When trilinear couplings are included we require that $a^{\rm SM}\lesssim 4\pi v/\sqrt{N}$ and $a\lesssim 4\pi m_{\phi}/N^{3/2}$, where $m_{\phi}$ is the typical mass scale of the scalars and $v\approx 174$ GeV is the vacuum expectation value of the SM Higgs.

Randomness

We introduce phenomenological randomness into the model through the mass matrix.333Randomness in the mass matrix was also studied in [26]. We use fully-populated matrices where each entry is drawn randomly from a Gaussian distribution. For $N$ real scalars the mass matrix is then symmetrized. After diagonalization such a random matrix has a known distribution for the mass eigenvalues in the large $N$ limit. In Eq. (1) we write the model after diagonalization.

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 $\mathcal{O}(1/\sqrt{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_{\phi}$. The modification of the couplings of the Higgs is governed by

[TABLE]

A weighted average of the CMS and ATLAS global signal strengths from 7, 8, and 13 TeV data [36, 37, 38] sets a limit of $\sin\theta<0.13$. The values of $a^{\rm SM}$ that we consider in Sect. 5 are all largely consistent with this bound.

Electroweak precision measurements receive corrections at one-loop [39]

[TABLE]

and are constrained at the permille level making them subdominant [40]. Above we have again imagined all scalars to have the same mass $m_{\phi}$.

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 searches, once the bound from Eq. (4) is taken into account, have no residual exclusion power [41, 42, 43, 44, 45]. However, the full 13 TeV dataset has not yet been analyzed and searches for resonant 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 $\lambda^{\rm SM}\sim\mathcal{O}(1)$ [46, 47, 48, 49].

4 Disordered Model Building

In this section we discuss the spectrum and eigenvectors of disordered mass matrices. We focus on $N\times N$ matrices with entries randomly drawn from Gaussian distributions. Our starting point is the Gaussian Orthogonal Ensemble (GOE), defined by

[TABLE]

where we use $M_{ii}$ to denote the diagonal entries of the mass matrix $M$ and $M_{ij}$ for the off-diagonal entries. In Eq. (6) the entries are drawn from the standard Gaussian distribution $\mathcal{N}(0,1)$ with mean $\mu=0$ and standard deviation $\sigma=1$.444The standard deviation of the off-diagonal entries is $1/\sqrt{2}$ rather than $1$ because we imagine that the matrix is symmetrized via $(M+M^{T})/2$. The GOE contains real symmetric matrices which correspond to models of real scalars. For complex scalars the analogous ensemble is the Gaussian Unitary Ensemble (GUE) that contains complex Hermitian matrices. Most of our results apply to both cases (where there are differences we specify them explicitly).

The natural generalization of the GOE, with broader physical applications, that we consider is

[TABLE]

where the diagonal entries are now drawn from $\mathcal{N}(\mu_{d},\sigma_{d})$ and the off-diagonal entries from $\mathcal{N}(\mu_{o},\sigma_{o})$. When we refer to the parameters as $\mu$ and $\sigma$ we are setting $\mu=\mu_{d}=\mu_{o}$ and $\sigma=\sigma_{d}=\sigma_{o}$ and drawing all entries from $\mathcal{N}(\mu,\sigma)$.

4.1 Summary of Main Results

The main result of this section is that in the large $N$ limit the spectral density of matrices with Gaussian-random entries follows a universal distribution called the Wigner semicircle distribution [50, 51, 52], which once appropriately normalized reads

[TABLE]

For finite $N$ and including the dependence on the parameters used for the Gaussian distribution the spectral density is

[TABLE]

The parameter $\beta$ is called the Dyson index and specifies the ensemble used. For the GOE we have $\beta=1$ and for the GUE we have $\beta=2$. The parameter $\mu$ does not appear in Eq. (9) as the distribution does not depend on $\mu$ with the exception of a single large eigenvalue at $\sim N\mu$. The spectral edges are at $\pm\sqrt{2\beta N}\sigma$.

In Fig. 3 we plot the eigenvalues of matrices with $N=20$ and $N=100$ along with the Wigner semicircle distribution.

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 $\sim 1/\sqrt{N}$. This is indeed the case as the sum of the squares of the components of the GOE eigenvectors follow a beta distribution with mean $\sim 1/N$ and variance $\sim 1/N^{2}$ [53].

We can infer several consequences from the Wigner semicircle distribution:

• •

All eigenvalues (except the largest one) are roughly contained in the interval $(-\sqrt{N}\sigma,\sqrt{N}\sigma)$. In particular if $\sigma\ll\mu$ they can be much smaller than the scale $\mu$ of the matrix entries.

• •

The typical eigenvalue spacing is $\sim\sigma/\sqrt{N}$.

• •

The lightest eigenvalue is $\sim\sigma/\sqrt{N}$.

The fact that $N-1$ of the eigenvalues fall in $(-\sqrt{N}\sigma,\sqrt{N}\sigma)$ even when $\mu\gg\sigma$ appears useful for model building. If one found a model where $\mu\sim M_{\rm Pl}$ while $\sigma\sim v$ it would appear that a hierarchy has been generated. However we have only rewritten $N-1$ Goldstone bosons in an unusual basis.

To see this, consider $N$ real scalars contained in a vector $\Phi$ with an $O(N)$-symmetric potential

[TABLE]

If we expand $V(\Phi)$ around the true minimum and take the VEV in the direction $\langle\Phi\rangle=(v,\ldots,v)^{T}$ the mass matrix of the physical degrees of freedom is precisely what we obtain in our examples by sending the standard deviation $\sigma$ to zero, i.e. a matrix with all equal entries. This matrix has $N-1$ massless eigenvalues which are the Goldstone bosons of SO($N$)/SO($N-1$). The large eigenvalue with mass squared proportional to $N$ is the radial mode. This explains why the typical scale of the eigenvalues is $\sigma$ which is explicitly breaking the symmetry, instead of $\mu$ that preserves it.

The last useful result on the spectrum is that taking a different mean for the diagonal elements $\mu_{d}$, shifts the spectral density from zero by the amount $(\mu_{d}-\mu)$

[TABLE]

If viewed in terms of Goldstone bosons this is a consequence of the explicit breaking of the $O(N)$ symmetry.

Before turning to the derivation of these results it is worth commenting 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.

[TABLE]

where $\langle\cdot\rangle$ 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 $\sigma^{2}\sim 1/N$ for our scalars, just from computing the free energy from the partition function $Z$

[TABLE]

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 can be long range.

Catalan Numbers

There is another interesting structural property of our large $N$ sectors with Gaussian-random mass matrices. The set of numbers

[TABLE]

known as the Catalan numbers, impacts both the shape of the eigenvalue distribution and the length of our decay chains. The moments of the Wigner semicircle distribution are the Catalan numbers

[TABLE]

A derivation of this result can be found in App. B while an explanation that uses an interesting 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

[TABLE]

with $N_{\rm max}$ determined by phase space. For $n\to\infty$ we have that $C_{n}\sim 4^{n}/n^{3/2}$ which means that the $n$th term in the average goes like $\sim(4pq)^{n}/\sqrt{n}$. Assuming that $(4pq)\sim\mathcal{O}(1)$ the average will go like $\sim\sqrt{N_{\rm max}}$. $N_{\rm max}$ is set by phase space and is $\mathcal{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 $\sim\sqrt{N_{\rm max}}$ continues to apply for higher $n$-body decays.

4.2 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.

The first step is to go from the joint entry distribution $\rho[M]$, in Eq. (7), to the joint eigenvalue distribution $\hat{\rho}[m_{1},\ldots,m_{N}]$ via the change of basis

[TABLE]

where $M_{D}={\rm diag}(m_{1},\ldots,m_{N})$. The joint distribution becomes

[TABLE]

Note that $\rho[M_{D}]$ is not the joint eigenvalue distribution because the Jacobian $|J(M_{D})|$ remains. The joint eigenvalue distribution is rather $\hat{\rho}[m_{1},...,m_{N}]=\rho[M_{D}]|J(M_{D})|\int dU$.

The metric on the space of symmetric (or Hermitian) matrices that defines $DM$ is induced by the product $M_{1}\cdot M_{2}={\rm Tr}[M_{1}M_{2}]$, i.e.

[TABLE]

where $ds^{2}_{M}$ is the distance on the space of matrices. The other differentials are

[TABLE]

For simplicity we start with the case $\mu=0$. Eq. (18) anticipates that the Jacobian of the transformation $|J(M_{D})|$ and $\rho[M_{D}]$ depend only on the eigenvalues. The proofs are given below and in App. A. If we think about $\rho$ 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$.

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

[TABLE]

To determine $\beta$ we can use dimensional analysis555We use square brackets $[\cdot]$ to indicate the dimensions of a quantity.

[TABLE]

from which we find that $\beta=1$ for the GOE and $\beta=2$ for the GUE. Now we are left with evaluating $\rho[M_{D}]$. If we take $\mu=0$ in Eq. (7) it is easy to conclude that in the large $N$ limit

[TABLE]

Combining this with the result for the Jacobian we finally obtain the joint distribution of eigenvalues

[TABLE]

where $Z_{N,\beta}$ is a normalization factor.

From Eq. (24) is possible to derive the Wigner semicircle distribution, Eq. (8), by solving a path integral [54, 55]. We can either use Feynman diagrams, as done in App. C, or use a saddle point approximation [50, 52]. Even without going through the derivation we can understand most of the results in the previous section just by looking at Eq. (24).

For example we expect the largest positive and negative eigenvalues to be $\mathcal{O}(\pm\sqrt{\beta N}\sigma)$ just from expanding $\hat{\rho}[m_{1},\ldots,m_{N}]$ around the largest eigenvalue $m_{*}$

[TABLE]

and taking the derivative $d\hat{\rho}/dm_{*}$ to find the maximum.

Thus far we have assumed $\mu=0$ but it can be shown that the results are valid also for $\mu\neq 0$. Consider splitting $M$ into a zero-mean matrix $M^{\prime}$ and a constant matrix $A$

[TABLE]

If we apply a unitary transformation $U_{A}$ that diagonalizes $A$ we find that $U_{A}AU_{A}^{\dagger}={\rm diag}(0,\ldots,0,N\mu)$. At the same time $M^{\prime}$ is rotated to $M^{\prime\prime}=U_{A}M^{\prime}U_{A}^{\dagger}$. The asymptotic distribution of eigenvalues of $M^{\prime\prime}$ is still the Wigner semicircle distribution and even the finite $N$ joint distribution is the same. The reason is that $\rho[M]$ is invariant under unitary transformations of $M$ and the Jacobian is the same for $M$ and $M^{\prime\prime}$ since they are both Hermitian.666For a derivation of this result see App. A

A similar argument can be applied to the case with a different mean for the off-diagonal entries: $\mu_{d}\neq\mu_{o}$. We can subtract the matrix $B$ where $B$ contains the mean of each entry. Then we further split $B$ into a constant matrix $A$ and $(\mu_{d}-\mu_{o})$ times the identity

[TABLE]

By diagonalizing $A$, $M^{\prime}$ is rotated as before, while the term proportional to the identity, $(\mu_{d}-\mu_{o})\mathbf{1}_{N\times N}$, is unaffected. Therefore we still have the Wigner semicircle distribution, but shifted from zero by an amount $(\mu_{d}-\mu_{o})$. In addition we have the large eigenvalue which is also shifted to $m_{\rm max}=N(\mu_{o}-1)+\mu_{d}$.

This concludes our heuristic derivation of results on the spectrum of GOE and GUE matrices. We now turn to another phenomenological possibility: having mass matrices that mix only $p$ nearest neighbors.

4.3 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$ nearest neighbors

[TABLE]

This scenario can arise, for example, from localization in an extra dimension. The eigenvalues of these matrices are spread over a smaller range than those in the GOE or GUE. $\rho[M]$ is the same as before leaving intact the Gaussian measure in Eq. (24), 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 find

[TABLE]

From this we find that

[TABLE]

In this case the eigenvalues are spread over an interval that we expect to be $\mathcal{O}(\sqrt{p}\sigma_{o})$ 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 [56], with mass eigenstates spanning approximately $p$ flavor eigenstates. This fact has already found several applications in the context of BSM physics [57, 58] and cosmology [59, 60, 61].

We do not explore band matrices in detail 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.

5 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, $\lambda^{\rm SM}$. The scalars are pair produced through this coupling (Fig. 5(a)) and can decay either through the same Higgs portal coupling (Fig. 6(b)) or through the hidden sector quartic $\lambda$ (Fig. 6(a)). For simplicity we take $\lambda^{\rm SM}=1/N$ and $\lambda=1/N^{2}$ for all the scalars.777We do not need this scaling to maintain perturbativity in all the cases that we study. We only use it as a convenient benchmark. The results presented here are not strongly affected by this choice of couplings.

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. 5(a)) and single production (Fig. 5(b)). The scalars can decay either as before (Fig. 6(a) and Fig. 6(b)), or through the two-body scalar channel (Fig. 6(c)), or to a pair of SM particles via mixing with the Higgs (Fig. 6(d)). For the quartic couplings we use $\lambda^{\rm SM}=0.1/N$ and $\lambda=1/N^{2}$. For the trilinears we take $a=m_{\rm min}/N^{3/2}$, where $m_{\rm min}$ is the lightest scalar mass, and $a^{\rm 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}\times m_{\rm min}$ and $a^{\rm SM}=10^{-5}\times 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 $\lambda^{\rm SM}$ and summed over all singly produced scalars in Fig. 7 (right) as a function of the trilinear $a^{\rm SM}$. We use a center-of-mass energy of 14 TeV and an integrated luminosity of $3~{}{\rm ab}^{-1}$. Cross sections were calculated using Madgraph 5 [62] with gluon fusion implemented via the Higgs Effective Theory module [63]. 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_{\rm tot}$ per event. In Fig. 8 we show $N_{\rm 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_{\rm 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 spectrum, 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 $\phi_{N}$ and the states closest in mass to it. When there is a large gap between $\phi_{N}$ and the next state, events with $N_{\rm tot}=2$ are favored because $\phi_{N}\phi_{N}$ production dominates the overall rate. Spectra for which the $N_{\rm tot}=2$ bin is small have states that are near in mass to $\phi_{N}$ so that their production rate is comparable to that of $\phi_{N}$. Furthermore, 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 $\phi_{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 250 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 pairs of SM objects are still the main avenue for discovery. However when pair production 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 $\phi_{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.

Let us now turn to the last scenario, where pair production dominates but $\phi_{N}$ can decay to a pair of SM particles. Interestingly for $N\lesssim 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 $\phi_{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 the Wigner semicircle distribution). The $WW$ decay width turns on rapidly for $m_{W}\lesssim m_{\phi}\lesssim 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 cannot 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.

6 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, light 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$ Dirac fermions with a Yukawa interaction

[TABLE]

This choice is rather appealing from the point of view of coupling this sector to the SM since it gives us a symmetry argument to forbid the $lH\psi_{\alpha}$ vertex leaving us with the interaction

[TABLE]

which gives the long decay chains discussed in Sect. 5. If we also take the mass of $\phi$ in a range that makes three-body decays within the dark sector ($\psi_{\alpha}\to\psi_{\beta}\phi^{*}\to\psi_{\beta}\psi_{\gamma}\psi_{\delta}$) 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^{\dagger}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\times N$ matrix $M$, with independent and identically-distributed entries drawn from a Gaussian distribution with zero mean and unitary variance the asymptotic distribution for the eigenvalues of $(M^{\dagger}M)/N$ is a particular case of the Marčenko-Pastur density [64]

[TABLE]

which denotes an accumulation of eigenvalues around zero. A diagrammatic derivation of this result can be found in [65]. Note that we have used $m$ for the eigenvalues of $M^{\dagger}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 $\sigma$, the asymptotic distribution for the spectral density of $(M^{\dagger}M)/N$ is

[TABLE]

Also in this case turning on a non-zero mean $\mu$ equal for all the entries does not affect the spectral density except for the appearance of one large eigenvalue of $M$ that is $\mathcal{O}(N\mu)$. We leave a more detailed discussion to future work.

7 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 $\mathcal{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.

Acknowledgements

We would like to thank Y. Bai, T. Cohen, B. Dobrescu, H. D. Kim, and G. Villadoro for useful discussions. RTD is supported by the U.S. Department of Energy under grant number DE-AC02-76SF00515. ML acknowledges support from the Institute for Advanced Study and from the Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. We thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work. This work was partially supported by a grant from the Simons Foundation (341344, LA).

Appendix A The Vandermonde Determinant

As discussed in Sect. 4 the joint entry distribution for real symmetric $N\times N$ Gaussian-random matrices with variance $\sigma$ and mean equal to zero is

[TABLE]

The joint eigenvalue distribution for the eigenvalues $m_{i}$ of the matrix $M$ is

[TABLE]

where $Z_{N}$ is a normalization factor.

The factor $|J|=\prod_{i>j}|m_{i}-m_{j}|$ is the Jacobian of the transformation between matrix entries and its eigenvalues and eigenvectors. It is known as the Vandermonde determinant. In this appendix we derive it in two ways.

The first is through a simple change of basis. Starting from the joint entry distribution we find

[TABLE]

To compute $J$, recall the definition of the metric on the space of symmetric matrices

[TABLE]

where we have used $d(OO^{T})=0$ which means $dO^{T}=-O^{T}dOO^{T}$. This also shows that $O^{T}dO$ is an antisymmetric matrix. We now have a metric tensor and we can use the square root of its determinant to obtain $J$.

Note that when we write $\rho[M]DM$ we are integrating only over the $N(N+1)/2$ independent variables in $M$. Given the form of $\rho[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. (38) to compute $J$. Combining Eq. (37) and Eq. (38) we obtain

[TABLE]

It is not hard to generalize these steps to the case of complex Hermitian matrices.

The second way to derive $J$ is as a Fadeev-Popov determinant. The correlation functions of gauge-invariant operators $O_{i}(A)$ in a Yang-Mills theory

[TABLE]

have the same structure as the expectation values of quantities that depend only on the eigenvalues of $M$

[TABLE]

if we identify $\rho[M]$ with the action $e^{iS(A)}$. In the matrix case a gauge transformation is a change of basis. From Eq. (35) we see that an orthogonal transformation on $M$ does not change $\rho[M]$

[TABLE]

and leaves us in the space of symmetric matrices $(OMO^{T})^{T}=OMO^{T}$. It is then easy to conclude that also in Eq. (41) the action and integration measure are gauge-invariant.

Following the Fadeev-Popov procedure we first define $\Delta(M_{D})$ as

[TABLE]

where $\delta[M-M_{D}]$ is an $N(N+1)/2-$dimensional Dirac delta function for symmetric matrices $M$. For complex Hermitian matrices it would be $N^{2}-$dimensional.

From Eq. (43) we compute $\Delta[M_{D}]$ for an infinitesimal orthogonal transformation $O\approx 1_{N\times N}+\delta O$ which means $DU\approx d\delta O$. We find

[TABLE]

So we have found that $\Delta(M_{D})$ is precisely the Jacobian of the transformation. For a complex Hermitian $M$, $|J|=\prod_{i>j}|m_{i}-m_{j}|^{2}$ where the square derives from the double integration on the real and imaginary parts of the elements of the unitary transformation.

The last step to derive Eq. (36) consists in multiplying Eq. (41) by the identity in the form given in Eq. (43) and note that $\int DU\delta[M-M_{D}]$ being the inverse of $\Delta(M_{D})$ cancels the factor of $\prod_{i>j}|m_{i}-m_{j}|$ in the integration measure. The Jacobian is then restored by $\Delta(M_{D})$.

Appendix B Moments of Wigner’s Semicircle

Recall the Wigner semicircle distribution:

[TABLE]

The odd moments of this distribution vanish by symmetry. Here we show that the even moments, $\mu_{2n}$, for integer $n$, are the Catalan numbers, $C_{n}$, where

[TABLE]

We compute the integral in polar coordinates

[TABLE]

Using a reduction formula

[TABLE]

the integral in Eq. (47) can be simplified to an integral over a single power of sine

[TABLE]

Finally using the identity

[TABLE]

one finds that

[TABLE]

Appendix 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

[TABLE]

Consider a Hermitian matrix $M$.888We show the Feynman rules for the GUE because they are slightly simpler than for the GOE. It has an associated Green’s function

[TABLE]

where $m_{i}$ are the $N$ eigenvalues of $M$. We are interested in the average over many realizations of $M$

[TABLE]

The averaging merges the poles into a cut spanning the support of $\rho(m)$. From the Green’s function one can find the spectral density $\rho(m)$ via the identity

[TABLE]

We compute $G(z)$ using the Feynman diagram expansion developed in [66] and discussed pedagogically in [55]. It is convenient to first introduce $G^{i}_{j}(z)$ as

[TABLE]

We expand $G^{i}_{j}(z)$ to find

[TABLE]

Odd powers of $M$ vanish in the average. The numerator of each term, explicitly, is

[TABLE]

where $Z$ is a normalization factor. Eq. (58) 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 [29]. 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 $\langle(M^{2n})_{j}^{i}\rangle$. Closed loops correspond 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 [67]. Therefore

[TABLE]

With Eq. (57) this implies

[TABLE]

which was shown explicitly in App. B.