Rydberg Composites

We introduce the Rydberg Composite, a new class of Rydberg matter where a single Rydberg atom is interfaced with a dense environment of neutral ground state atoms. The properties of the Composite depend on both the Rydberg excitation, which provides the gross energetic and spatial scales, and on the distribution of ground state atoms within the volume of the Rydberg wave function, which sculpt the electronic states. The latter range from the"trilobites", for small numbers of scatterers, to delocalized and chaotic eigenstates for disordered scatterer arrays, culminating in the dense scatterer limit in symmetry-dominated wave functions which promise good control in future experiments. We characterize these scenarios with different theoretical methods, enabling us to obtain scaling behavior for the regular spectrum and measures of chaos and delocalization in the disordered regime. Thus, we obtain a systematic description of the Composite states. The 2D monolayer Composite possesses the richest spectrum with an intricate band structure in the limit of homogeneous scatterers.


I. INTRODUCTION
Ultra long-range molecules composed of a Rydberg atom and a ground state atom, colloquially known as trilobites, were proposed in 2000 [1]. Soon thereafter theoretical explorations regarding the possibility of polyatomic molecules involving several ground state atoms followed [2,3]. The experimental verification of Rydberg molecules in 2009 [4] also confirmed accidentally the existence of trimers [5]. Since then, interest in Rydberg excitations beyond isolated atoms has rapidly branched out into quite diverse scenarios. These include the replacement of the ground state atom in the original trilobite dimer by larger and more complex systems, e.g., one or more polar molecules [6][7][8][9], the (re-)discovery of Rydberg excitations in solid state systems [10], and a large variety of excitonic Rydberg dynamics in the gas phase [11][12][13], just to name a few. For the increasingly dense gases now achievable in experiments, one can elegantly describe this system as a Rydberg excitation dressed by ground state atoms from the gas. In fact, recent experiments exhibit spectral features corresponding to polyatomic molecules containing up to five ground state atoms [14][15][16], and mean-field shifts in the spectrum reveal this polaronic behavior involving the coupling of many hundreds of atoms to the Rydberg electron [17]. One may wonder how many ground state "scatterer" atoms within the volume occupied by the Rydberg wave function can a trilobite molecule tolerate. A recent study found that under certain conditions the formation of trilobites actually thrives in a dense gas, which is counter-intuitive at first glance [18].
What is lacking is a systematic approach which connects the trilobite regime with a few scatterers to the regime of very dense scatterers, although the phenomena just described suggest that Rydberg excitations immersed in dense and structured media might have very interesting properties. The present investigation opens a new venue for Rydberg Composite systems along this way, which involve many hundreds of atoms in a struc-tured environment coupled to a single Rydberg atom. These Composites can be formed, for example, by exciting an atom in a 1, 2, or 3-dimensional optical lattice to a Rydberg state which envelops many atoms on surrounding sites. We present a systematic and detailed investigation of this Rydberg Composite and provide its properties as a function of principal quantum number ν, lattice constant d, and fill factor F of lattice sites.
With the Rydberg composite we change the perspective from the molecular one -using chemical approaches to characterize polyatomic trilobites via Born-Oppenheimer potential surfaces, rovibrational couplings, etc [2,19] -to a condensed matter one, emphasizing generic scaling principles, gross structure, and properties associated with the high density of states obtained here. This allows us to approach systematically dense atomic environments. Indeed, we will see that towards the limit of homogeneous filling a band-like structure of the spectrum emerges. Moreover, the unique property of a Rydberg electron bound to an isolated atom with a singular point of infinite density of states (DoS) at the ionization threshold lim ν→∞ E ν ≡ −1/(2ν 2 ) = 0 and full degeneracy makes such a Rydberg Composite an interesting object to study, as the distribution of scatterers can break the degeneracy in a controlled, yet flexible, way. We will identify non-trivial scaling properties as a function of ν. They allow us to connect the situation at finite ν with threshold ν → ∞. Finally, the Composite's key properties are derived analytically in the homogeneous limit, while random matrix theory is used for the irregular part of the spectrum.
We will also explain how a planar environment breaks the symmetry of the Rydberg Composite and leads to much richer spectral structures as compared to a wirelike (one-dimensional) or crystal-like (three-dimensional) atomic environment. Hence, we put emphasis on a planar sheet of atoms arranged in a lattice containing a Rydberg excitation as an exemplary Rydberg Composite whose experimental realization is facilitated by the routine creation of two-dimensional optical lattices [20]. This paper is structured as follows. Sec. II provides the theoretical background. In II A we introduce a generic Hamiltonian which represents a broad class of systems consisting of an excited object coupled to localized scatterers. Sec. II B defines our specific realization of this Hamiltonian: the Rydberg Composite in D dimensions. Sec. II C details the scaling properties of this Composite system in one, two, and three dimensions. In Sec. III we introduce the phenomenology of the Composite in the three different lattice geometries, investigating both the DoS (Sec. III A) and exemplary wave functions (Sec. III B). In section IV we focus on the homogeneous density regime where the system can be studied analytically to obtain a clear intuitive picture of the system, its band structure, and the resulting scaling laws. Section V investigates the inhomogeneous regime, using statistical measures derived from random matrix theory to reveal that it exhibits quantum chaos. Sec. VI discusses potential experimental realizations, and Sec. VII concludes with further perspectives and implications. Throughout we adopt atomic units.

A. Generic Hamiltonian
We begin with a generic description of our system, which is composed of an electron with position r and momentum p in the presence of a central potential V (r) and a collection of point-like scattering objects with positions following a distribution of particles ρ( x). This scatterer arrangement can correspond to either a structured geometry or a disordered environment, i.e. that found naturally in an optical lattice or in an ultracold gas, respectively. Although the electron wave function is fully three-dimensional (3D), the dimensionality of the scatterer geometry can be lower, for example as in a onedimensional (1D) chain or a (2D) random gas. The scatterers interact with the electron via the potential U ( x, r), which destroys the spherical symmetry of the central potential V (r) and, in general, makes the system classically chaotic. We assume a frozen-gas scenario, consistent with the ultracold temperatures of such a system, and neglect the motion of these scatterers. The electronic Hamiltonian is therefore This generic Hamiltonian has been studied in several contexts over the past decades, with examples ranging from two-dimensional (2D) quantum dots [21,22], quantum billiards [23], Coulomb systems [24,25], perturbed harmonic oscillators [26], and Bose-Einstein condensates in a dimple potential [27], to name just a few. The electronic wave function for vanishing U separates in spherical coordinates: is a spherical harmonic. Therefore, mutual eigenstates of the angular momentum operator and the Hamiltonian satisfy L 2 |νlm = l(l + 1)|νlm , L z |νlm = m|νlm , and H|νlm = E νl |νlm . Any central potential possesses azimuthal symmetry and hence has 2l − 1 degenerate |νl states. In the next section we consider the Coloumb potential, V (r) = −1/r, which has an additional symmetry: it conserves the Runge-Lenz vector A = p × L −r, leading to a particularly large degenerate Hilbert space in each manifold ν. Scatterers will lift this degeneracy. Special scatterer geometries, however, may be able to restore this degeneracy in the Rydberg Composite.

B. Implementation for the Rydberg Composite
The Rydberg atom is a major workhorse of modern atomic physics; here, when embedded in an ultracold medium of neutral atoms, it provides an ideal physical realization of the Hamiltonian (Eq. 1). For an alkali atom, this means that V (r) is a Coulomb potential for r larger than a few atomic units. The deviation at small r, typically set by an empirical model potential, includes the interactions with the other atomic electrons. This leads to energies E νl that are non-degenerate for different l values and noninteger values of ν, the principle quantum number. However, as l increases, the wave function's overlap with this short range region decays rapidly and E νl → − 1 2ν 2 , where ν is an integer. Typically only the states with the three or four lowest l values deviate appreciably from the hydrogenic Rydberg spectrum. The overwhelming majority of states behave as in hydrogen, and therefore for simplicity we consider only a hydrogenic spectrum here with ν an integer. For the interaction between the surrounding ultracold atoms -the localized scatterers -and the electron we use the Fermi pseudopotential [1,28], which is straightforward to implement and manipulate. The strength of each scatterer's contribution is given by the energy dependent s-wave electron-atom scattering length a s [k x ] [29]. This simple contact potential is a reasonable approximation since a neutral atom in its ground state is a highly localized and isotropic perturbation when compared to the Coulomb potential and Rydberg wavelength. It therefore imparts only an s-wave phase shift onto the Rydberg wave function via elastic scattering, as characterized by the scattering length. The energy dependence of this process is set by the semiclassical electron momentum, k 2 Schematics of the three scenarios we consider: (a) a linear chain (1D), (b) a monolayer (2D), and (c) a cubic lattice (3D). In each panel the black spheres represent scatterers sitting on lattice sites and the red lines give the lattice spacing, d. The missing scatterers in (b) represent a situation with non-unity filling. The volume of scatterers situated within the Rydberg wave function is represented by the blue circle in (b) and by the sphere in (c). The exemplary densities shown give different representations of the Composite's electronic wave function in these three scenarios. In (a) the strongly perturbed wave function is shown with three surfaces of constant density, revealing the exotic nature of these wave functions. The full 3D contour is cutaway in front to reveal the interior structure. Fig. 3 provides details of the plot parameters. Panel (b) shows a contour of the angular dependence of a typical circular state, which plays a crucial role in the 2D-Composite's properties (see Sec. IV). Panel (c) shows a cartoon of the Rydberg atom, illustrating that the spatially-varying probability cloud spans many lattice sites of the 3D lattice.
The eigenvalues E i and eigenstates |Ψ i = νlm c (i) νlm |νlm of this Hamiltonian are parameterized by the distribution of scatterer locations, ρ( x). We are interested in scatterers in a lattice configuration, and hence choose where the scatterer positions R i = d D j=1 n ijêj are located at lattice positions described by unit vectorsê j , the lattice spacing d, and a set of D × N D integers n ij . By excluding some values i we can implement partial filling, defined by the fill factor F . Eq. 3a relies on two more approximations: the scattering length is energy-independent and the basis is truncated to only a single ν-manifold. We demonstrate in Sec. IV that these approximations are increasingly accurate at high ν. They have only minor quantitative effects on the main conclusions of our study but allow us to obtain analytical formulas and clear scaling behavior. Since we consider only a single ν-manifold, in the following discussion we will set − 1 2ν 2 to zero. Hence, the spectrum of Eq. 3a with the distribution from Eq. 3b is obtained by diagonalizing

C. Scatterer induced properties in D-dimensional lattices
The properties of the Rydberg Composite, defined by the Hamiltonian in Eq. 3a, depend on the properties of both the unperturbed electronic states, |νlm , and the scatterer distribution, ρ( x). The Rydberg atom's size, density of states, and wavelength are determined by its principal quantum number ν, while ρ( x) depends on the desired lattice geometry, lattice spacing, and filling realization. In this section we delineate the important quantities for one, two, and three-dimensional scatterer configurations.
Although the spatial scale of the lattice greatly exceeds the size of the Rydberg wave function, not all scatterers perturb the Rydberg states since the electron-atom interaction is highly localized. Its strength, within the Fermi approximation, is determined by the electronic density directly at the scatterer position. The Rydberg volume is finite with a radius r 0 (l) ≈ a l ν 2 , where a l decreases from a l ≈ 2 for the l = 0 state to a l ≈ 1 for the l = ν − 1 states. Numerically, we consider scatterers inside a radius r = aν 2 with a > 2 to guarantee that all wave function amplitudes are exponentially small, and hence contribute no energy shift, at the boundary of this volume. The number N D of relevant scatterers is then determined by the volume V D of the intersection of the Rydberg wave function and the lattice. In this way, even for an infinite lattice, we can truncate its effect to only that caused by the N D individual scattering potentials in Eq. 3a.
A 1D array of scatterers corresponds to a chain lattice, see Fig. 1a. The relevant 1D volume is V 1 = 2aν 2 , and N 1 = 2aν 2 /d scatterers lie within this volume. The corresponding volume for a 2D lattice is the area of the projec-tion of the Rydberg volume into the plane, V 2 = πa 2 ν 4 , and hence the number of scatterers is N 2 = πa 2 ν 4 d 2 . In 3D we consider a cubic lattice of scatterers, and so the relevant volume is the entire Rydberg volume, V 3 = 4 3 πa 3 ν 6 , containing N 3 = 4πa 3 ν 6 3d 3 scatterers. The N D values given here are valid only in the N D 1 limit where edge effects due to the incommensurate spherical and cartesian geometries are negligible.
The scatterer configuration also influences how many of the degenerate states of the Rydberg manifold are shifted. As a general rule, each scatterer splits away one state from the degenerate manifold until the geometry induced limit B D is reached. In a generic 3D scatterer array this limit is given by all states of the manifold, B 3 = ν 2 , while B D < ν 2 in 1D and 2D. To determine B D for each case we select a convenient quantization axis and identify the Rydberg states not affected by the delta-function potential (Eq. 2). In 1D, we set the quantization axis parallel to the linear chain of scatterers. When r → Rẑ, most angular wave functions vanish on the quantization axis since Only m = 0 states experience a shift, and hence B 1 = ν is the total number of m = 0 states. For the 2D case we set the quantization axis normal to the plane and evaluate the angular wave functions at θ = π/2. The Legendre polynomials with argument cos(π/2) = 0 are The plane is transparent to the Rydberg states possessing a node in the plane. Therefore, With the help of B D we can define a third quantity, the characteristic lattice spacing d D such that N D ≈ B D for that geometry. This spacing heralds the onset of the density shift regime where additional scatterers cannot split away new states since the ν-manifold is saturated. They instead contribute linearly to a mean-field energy shift, consistent with the conclusion drawn from the original applications of Fermi's pseudopotential [30,31]: the mean-field effect of the interaction of the Rydberg electron with the scatterers is an energy shift proportional to the electron-atom scattering length and to the scatterer density. The values for this characteristic length, along with the other values V D , N D , and B D , are given in table I. From these characteristic properties we can assess the behavior and crude scaling with ν of the Rydberg Composite for a given scatterer geometry. Notice that for D ≤ 2 the critical lattice spacing is linear in ν, but follows ν 4/3 for the 3D case.  This analysis suggests that for sufficiently large number of scatterers N D we will obtain B D non-zero eigenvalues upon diagonalizing H within a ν-manifold, and as a function of decreasing d these eigenvalues will grow (on average) linearly with the number of scatterers. In order to remove this asymptotic shift we normalize the total energy shift by N D . Furthermore we measure energies in units of (2π|a s |) −1 in order to remove the numerical prefactor from the potential matrix (Eq. 4), and hence eliminate the material-dependent value of the scattering length from our calculated energy shifts. Finally, since N D depends on the arbitrary (provided it is sufficiently large) choice of a, we scale the energy shifts of the Ddimensional lattice by a D to eliminate this scale choice. Of course, in the limit a → ∞, this choice removes all dependence on a and we can report scaled energiesẼ defined in terms of the un-scaled eigenvalues E viã whereṼ D is the volume of a D-dimensional sphere with unity radius. We now investigate the behavior and properties of the Rydberg Composite by computing its spectrum for each geometry.

III. PHENOMONOLOGY OF THE RYDBERG COMPOSITE
The spectrum of a Rydberg atom immersed in a structured neutral medium depends both on the Rydberg principle quantum number and on the different realizations of the lattice. We parameterize the latter by its filling factor F , the percentage of filled lattice sites, and by the lattice spacing d. We focus first on unity filling factor so that we can introduce the essential quantities useful in characterizing the Composite's properties. In section V we will remove this restriction and study fractional filling.
We first study the density of states (DoS). It reveals more about the global spectral properties than individual energy levels, and provides a useful guide to regions of interest to focus on in finer detail. In a second step, guided by the features seen in these DoS, we study the wave functions corresponding to various paradigmatic states. The structure present in these wave functions provides additional investigative tools to understand the spectra. Since the 2D monolayer leads to the richest structure in the dense lattice limit, we focus on that geometry.

A. Density of States
We show DoS in Fig. 2 for the lattice geometries depicted in Fig. 1. We observe that all B D eigenstates converge to constant limits for d → 0, as anticipated. Intriguingly, we find that the asymptotic value differs remarkably across the three geometries. For the 1D and 3D scatterer geometries the shifted eigenenergies become degenerate again as N D → ∞, albeit at a large overall energy shift relative to the zero-scatterer degenerate manifold. In contrast, eigenenergies in the 2D geometry remain non-degenerate even in the infinite density limit, instead developing three main features (see also the spectrum Fig. 5): a nearly continuous and quasi-uniformly spaced distribution of energy levels within a few "bands", the formation of a large "band gap" that persists even up to relatively large lattice constants, and the formation of a large peak in the DoS in the upper part of the spectrum.
In all three lattice geometries, as d increases the DoS becomes challenging to interpret due to the increasing number of non-degenerate energies. In general, the spectrum diffuses. In 1D, the degenerate band is depleted as individual states split discretely away with increasing d. This process does not occur symmetrically with respect to the degenerate band. In 3D, all states begin to split apart at approximately the same value of d and the perturbed band dissipates far more rapidly than in 1D; this process also occurs symmetrically about the homogeneous energy asymptote. In 2D the states are not degenerate in the d → 0 limit. For increasing d, states higher in the energy band begin to disperse linearly in d, revealing a clear energy dependent transition between the indistinguishable (d ≈ 0) and distinguishable scatterer case. In all three geometries, oscillations in the energy levels mimic the oscillatory nature of the Rydberg wave function, which is imposed quite directly onto the energy levels via the contact potential. The "spaghetti" nature of the energy levels in the large d regime reveals the presence of both, real and avoided level crossings if d is taken as an adiabatic parameter. Real crossings are possible since the electronic states mirror the lattice symmetry, and therefore can be grouped according to the irreducible representations of the nuclear point group for that lattice. We have confirmed that, in the 2D lattice case, the DoS can be computed independently for each of the five irreducible representations of the C 4v point group, following the description of Refs. [2,19]. This is discussed in more detail in Appendix C. Finally, as d grows further, the DoS (not pictured) collapses gradually back into a highly degenerate peak at zero energy as the number of scatterers falls below B D .

B. Wave function characteristics
We now present a representative sample of the wave functions giving rise to these DoS in 1D and 2D, which are particularly amenable to this treatment since all relevant information can be gleaned and easily visualized with three-dimensional contour plots (1D) or the z = 0 slice through the electron density (2D). From these wave functions we begin to see the underlying structure of the perturbed system and how it might lead to the emergence of the structured, non-degenerate bands in 2D band structure rather than a single, fully degenerate band in 1D and 3D. Although our focus now is descriptive, merely commenting on the appearance and classification of these wave functions, we will use these observations in the following section to develop quantitatively accurate approximations which lead to a full interpretation of the Rydberg Composite's properties.
3. Electronic densities of a 1D chain of scatterers for ν = 30. The yellow, green, and blue surfaces correspond to contours at wave function densities spanning factors of 10. The full three-dimensional surfaces are cut away in front to reveal the interior structure. The left column gives the density for the most deeply perturbed state, while the right column gives the density just above the degenerate band.

1D lattice wave functions
We present, for four different lattice spacings, two representative wave functions for the ν = 30 1D Rydberg Composite. On the left we show the state with the largest energy shift, while on the right we choose a state slightly higher in energy than the degenerate band limit, i.e. one of the states visible in Fig. 2a just above the middle band. These wave functions visually forge the connection between Rydberg Composites and "trilobite" molecules [2,3,19,32]. At large d, shown in the bottom row, the wave function is a mixture of many l states, leading to strong localization on scatterer positions. In scenarios such as this where scatterers are separated by distances greatly exceeding d D , the wave function tends to localize on only a subset of the scatterers and effectively ignore the rest. In this way it maximizes the overlap between the Rydberg electron and the lattice, and ensures orthogonal wave functions. As d decreases these states eventu-ally begin to resemble the hydrogenic basis states and localize less severely on a symmetry-imposed collection of scatterers, as the Rydberg wave function increasingly cannot distinguish scatterers lying closer together than its wave length. Unfortunately, these wave functions do not as yet reveal with any clarity why the infinite density limit of this 1D-Composite is again an energetically degenerate system. A key reason for this uncertainty is, in fact, their degeneracy: degenerate eigenstates obtained via a numerical diagonalization will in general be arbitrary superpositions of the, linearly independent, states. It is thus impossible to identify any possible good quantum numbers or selection rules from these wave functions without investigating some other observable. In principle, this could be done by applying a magnetic field to break apart the degeneracy at large scatterer density. For our present purposes we can turn instead to the 2D-Composite, which is fundamentally non-degenerate in this limit and may reveal through its wave functions the underlying structure of the 1D case.

2D lattice wave functions
Since only the electronic density in the z = 0 plane contributes to the energy shifts, it suffices to examine |Ψ(x, y, 0)| for the 2D-Composite. We first consider 2D-Composite wave functions with a fully filled 2D lattice and vary the lattice constant. In Fig. 4(Top) we show the wave functions corresponding to the first three oddnumbered eigenenergies starting from the lowest one. For large enough d the electron density obeys one of the discrete symmetries permissible by the C 4v point group, and partially localizes on only a subset of the available scatterers. The behavior of this localization and its effect on the energy level structure likely warrants future study. As d shrinks further, the electron density evolves into a distinctly circular shape. By the lowest d shown (d = 20), these three eigenstates have seemingly converged into "circular" states. By a strict definition, a circular Rydberg state has l = m = ν − 1; here we employ a broader definition meaning a state with high l and |m|, but with only a small difference l − |m|. The second and forth eigenstates identically resemble the first and third, respectively, showing that the ±m states are equivalent and degenerate in this limit To confirm that these eigenstates do not arise due to some coincidence in the symmetry-adapted wave functions or fortuitous overlap with the lattice grid, we next consider a lattice with a small lattice constant d d D but with varying fill factor F . At extremely low F (first column of Fig. 4(bottom)), having only a very few scatterers, we see that the few non-degenerate eigenstates are basically independent trilobite dimers between the Rydberg core and each individual scatterer. As F increases the number of scatterers increases rapidly and the wave function becomes rather chaotic in appearance, exhibiting no clear structure. In some instances it local- izes asymmetrically about statistical fluctuations in the random scatterer distribution where small clusters form spontaneously.
As before, when the series progresses towards complete filling, the density resembles more and more a circular state, thus confirming that the appearance of such states depends more on the total density of scatterers relative to the number fluctuations caused by random fill factors than on the underlying lattice symmetry. Once fluctuations and correlations in the scatterer density are unresolved by the Rydberg wave function, any choice of random fill factor is essentially indistinguishable and the result from the F = 1 case is reached.
C. Role of wave function character on the 2D spectrum Both ways of increasing the scatterer density described above lead to the following conclusion in the high density limit: the wave functions become increasingly circular in character, implying that they become approximate eigenstates ofL z . This is to be expected in the limit of a totally homogeneous lattice, where H commutes withL z due to the cylindrical symmetry. Of greater interest is the fact that the energies of these states are also apparently sorted by the level of circularity, as states with the most circular character fall to the bottom of the energy bands. A useful diagnostic to analyze the evolution of the wave functions as d or F changes is the participation ratio of m states, which ranges from 1 for a state proportional to δ mm to 1/ν for a state mixed uniformly among m sub-levels. In Fig. 5 we show the 2D-Composite's eigenspectrum as a function of d. Although this conveys very similar information as the DoS plot in Fig. 2, coloring the eigenstates by P m reveals additional structure in these energy levels that can be linked to the wave function. The P m distinguishes many self-similar and repeating substructures that were not evident in Fig. 2b. Several "bands" of states with a similar functional dependence on d and pattern of P m are visible, separated by the large energy gap that was clear in the DoS as well. These bands become indistinguishable towards high energy and converge into the region of high degeneracy seen in Fig. 2. The clear transition between states with P m ≈ 1, which have m as a good quantum number, and those which are strongly mixed helps differentiate these bands even when they start to overlap. This transition is well predicted by a critical lattice spacing (d c ) defined in Sec. V and shown for the first three bands in Fig. 5 as black curves. Figure 5 shows that the trend towards circular states in the few eigenstates presented in Fig. 4 is emblematic of a more general behavior: in the d/ν 1 limit H commutes withL z and hence all wave functions have m as a good quantum number, of which the circular states are a small subset. From the wave functions in Fig. 4 we see also that the lowest states of the energy bands approximately conserve l as well, with maximal or nearly maximal values of both quantum numbers.
We will devote much of the reminder of this paper to the 2D-Composite, since embedding the Rydberg excitation in a planar environment constitutes a new scheme in ultracold Rydberg physics with a rich and intricate behavior. In particular, we will explicate the physics underlying these (so far phenomenological) observations: the formation of energy bands separated by a single dominant band gap and the relatively simple character of the underlying wave functions. We elucidate the link between this structure and the structure underlying the Hamiltonian matrix. The mathematical tools used for this task also enable us to understand why all states become degenerate in a 1D and 3D homogeneous environment but not in a planar one. Furthermore, these tools prove useful as a launching point for our later investigation of disorder in dilute lattices with random filling and lattices with larger d.

IV. RYDBERG COMPOSITE PROPERTIES IN THE HOMOGENEOUS DENSITY LIMIT
Our investigation of the 2D-Composite spectrum starts with the observation that, below a certain lattice constant d < d c , the Rydberg wave function can no longer resolve individual scatterers. The lattice then appears homogeneous, and the phenomenology of previous sections has shown that the spectrum becomes constant. In section V we will clarify this coarse-graining concept further for the disordered scenario. In the present section, we will take it as fact that this coarse-graining is physically relevant and use it to approximate the discrete lattice of scatterers with a continuous plane of homogeneous density. In this way we characterize the system's properties for the d ν region of Figs. 2 and 5, which is then crucial to properly situate our analysis for intermediate cases with d > d c or F < 1.
A. 2D Monolayer: emergence of a band structure The replacement of the discrete lattice with a homogeneous distribution coincides mathematically with the replacement of the summation in Eq. 4 with an integral. In the scaled energy units this replacement must include also a factor V −1 2 , and the matrix elements become with integration over the entire plane. Using the spherical coordinate representation of these wave functions, the three contributions to the matrix elements are given by the product of an integral over ϕ, a radial overlap integral, and the projection of the spherical harmonics into the plane, P lm,l m = N lm P m l (0)N l m P m l (0), where and P m l (cos θ) was given in Eq. 6. As expected, integration over ϕ imposes a block-diagonal structure in m on this matrix, since the homogeneous scatterer limit is isotropic. Eq. 10 therefore yields of upper bands begin to overlap the flat low-|m| regions, and we find in this overlapping region that each band begins along the essentially continuous lineẼ = − m √ 2πν 3 . To understand the formation of these bands as well as to obtain analytic results for the eigenspectrum, we study the matrix elements of each m-level block of the Hamiltonian. As one might surmise from the wave function study in the previous section, these sub-blocks appear to be diagonal-dominated for moderate to high l and m. This can stem from two influences. First, the off-diagonal couplings tend to be around one order of magnitude smaller than the diagonal matrix elements. This is because, although the radial overlap integral does not have a rapid dependence on l, the u νl−2 (R) wave function has a node nearly at the maximum antinode of the u νl (R) function, and hence the integrand of R (1) νl,νl is small compared to R (1) νl,νl . An additional and more critical contribution stems from the spherical harmonic projections onto the plane, which are illustrated with surface contour plots in Fig. 7 for these circular states. Each column of this figure contains all the states within a single m-block, starting with m = ν − 1 = l on the left and decreasing by 2 (because of the selection rule of Eq. 6) with each step to the right. The orbital angular momentum l decreases by 2 with each vertical step up from its maximal value, l = ν − 1, in the bottom row. Since l cannot be less than |m| the size of each sub block increases with decreasing m. There is no coupling between columns, which come from different m-blocks. Clearly, the states within each column have dramatically different overlap with the z = 0 plane since each drop in l pushes an additional lobe out of the plane. These contribute nothing to the total energy shift and are essentially "wasted" probability. In contrast, along the diagonals marked by red the wave functions are nearly identical. Moving up the diagonal swaps an angular lobe into a (not pictured) radial lobe, which has a negligible impact on the overlap with the plane compared to the loss of pushing an entire set of lobes out of the plane. Within each m-block, therefore, the diagonal elements will have large energy separations, and hence effectively decouple. The resulting states that share similar qualities are spread over many blocks and correspond to the series along the diagonal in Fig. 7. We label elements in these series with the values b = 1 . . . ν − 1 and k = 1 . . . ν − b. These numbers label the diagonals and the state in a given diagonal, respectively. Note that from the construction in Fig. 7 we have l = 2b − 2 + m and k = ν − l. The Fig. 7 inset displays the overlap P lm,lm for various b as a function of l. Each b-band smoothly changes with l, but the difference between b-bands is large, especially so for the lowest b values.
The appearance of bands in the eigenspectrum (see Fig  6) is essentially a consequence of the approximate b-bands of similar states described in Fig 7. The identity b = β holds exactly when the m-block matrices are diagonal, i.e. in the asymptotic wings of these bands where the dispersion becomes approximately linear (see Fig. 6) and where l ≈ m ≈ ν. Since in this limit the eigenenergies are obtained analytically and the real bands β coincide with the approximate bands b, we focus now on the behavior of these linear wings. If we consider only the diagonal matrix elements of Eq. 14, we obtain the energy levels for the ∼ 10 lowest energy levels in each band of states to a few percent accuracy, confirming that the diagonal approximation is appropriate. More importantly, by switching to the band numbers b and k we can gain further intuition into the true bands, labeled by β, seen in Fig. 6. In the high ν limit we obtain the linear dispersion relatioñ In particular, band b ≈ β starts at the energỹ For b ν the energies scale as ν −11/2 . In particular, the lowest energy lies atẼ = ν −11/2 π −3/2 . At higher b we use the limiting form of the Γ functions, to obtain The band's lower edges become more closely spaced in energy due to this 1/ √ b dependence. The level spacing within a band is given approximately by Taking this width as approximately constant over an entire band and taking the number of states per band to be ∼ ν, we find that the width of each band is approximately On the other hand, the spacing between band minima is approximately Within this crude series of approximations we find that Due to this decreasing gap between bands relative to their own widths the bands begin to overlap, leading to the region of high energy density seen in Figs. 2b and 5, and apparent in Fig. 6. As the bands overlap with increasing b, the expression for the band minimum Eq. 17 tends towards ν −6 , as b ∼ ν. Apparently, as the energy-level structure transitions from the "band" type into this denser structure of many overlapping bands, the functional form of the energy scaling changes from ν −11/2 to ν −12/2 . B. Scaling laws for scatterers in D-dimensions

2D scatterers
The preceding analysis showed that the energy spectrum of the 2D-Composite exhibits two different scaling behaviors as a function of band number. The energies in the lower bands scale as ν −5.5 , but they scale as ν −6 in the upper bands due to a mixture of overlapping bands and deviations from the diagonal approximation for small values of l and m.

1D scatterers
Applying this same analysis to our 1D and 3D configurations leads quickly to the results that all states experience an identical energy shift. In 1D, the expression equivalent to Eq. 10 is Curiously, this radial matrix element vanishes when l = l [33]: This leads to a simple expression for the energy shifts, The Rydberg levels are identically affected by the scatterers and remain degenerate. Explicitly inserting the volume, we haveẼ The scaled energies scale as ν −5 , decreasing slower with increasing ν than the 2D-Composite energies. (a)

3D scatterers
For a homogeneous structure in 3D the matrix elements are even simpler, This is the normalization integral, and thus all Rydberg states are again degenerate, but with a global shift, These scale as ν −6 , more slowly than the 2D-Composite. The Rydberg Composite spectrum thus obeys a powerlaw scaling behavior ν f , where f = −5 in 1D, f = −11/2 in 2D, and f = −6 in 3D. Despite the fact that the monolayer gives an energy scaling intermediate between the other geometries, it leads to a non-degenerate, highly structured dense limit which is totally distinct from the 1D and 3D-Composites.

Interpolation of low and high band edge scaling for 2D scatterers
We now explore the DoS scaling in the 2D case in further detail to arrive at a universal DoS for 2D Rydberg Composites in the homogeneous limit. We compute a smooth density of states, where F is a convolution function for the discrete data, i.e. a Gaussian or a box function centered at E i and having width σ. We focus first on the "band" region, where the eigenfunctions are to a good approximation labeled by integers b and k, and which scale as ν −11/2 . Since the number of states in the bands increases approximately linearly with ν, we rescale the widths also so that they decrease linearly in ν. Fig. 8a shows the resulting DoS for three ν values using a Gaussian distribution for F . The agreement between different ν is excellent in this band region, breaking down as energy increases. In grayscale we overlay each band separately, showing how the total DoS is built up from these, and in particular how the overlapping bands create the saturation point in the DoS and the eventual onset of the ν −6 scaling law. Fig. 8b shows a DoS with the ν −6 scaling, appropriate to the "overlap region" where the diagonal approximation breaks down.
The widths now decrease as σ/ √ ν to obtain a smooth function. Some technical details involved in these figures are discussed in Appendix D. Fig. 8c presents a scaling that smoothly interpolates between these two regimes as a function of E. It allows us to construct a "universal" density of states for the 2D-Composite, independent of ν. Details of this process, which uses a hyperbolic tangent to map the relevant scale factors, widths, and normalizations between these two regions as a function of E, are provided in Appendix D. The DoS shown in Fig. 8c confirm that this scaling is indeed universal, as the DoS for the three different ν levels are essentially indistinguishable.

V. EVOLUTION OF THE SPECTRUM FOR DECREASING DENSITY OF SCATTERERS
We have thus investigated the nature of the spectrum in the limit where the lattice cannot be resolved by the Rydberg wave function. For the 2D-Composite this led to a non-trivial spectral density with band-like structures, scaling laws, and Rydberg wave functions quite different from the ones known from atoms or molecules. In this section we study the characteristics of this system at lattice spacings large enough to be resolved by the Rydberg wave function. To this end we define a threshold lattice spacing, d c , below which the Rydberg wave function can no longer resolve the scatterers and which therefore replicates the homogeneous limit of scatterers, formally only reached for d = 0.

A. Transition to homogeneous density of scatterers
Several circumstances complicate a rigorous definition of d c . First, the electron's wavelength varies spatially: in the radial direction it increases quadratically, while in the angular degree of freedom it is strongly l-dependent. Secondly, since the potential depends non-linearly on the wave function amplitude at the locations of the scatterers, it is not clear from the onset at what length scales scatterers can be resolved.
We have already seen that the homogeneous limit in the 2D case is heralded by wave functions which are diagonal in m. Near the bottom of each band these states are also approximately diagonal in l. Such a wave function has 2m angular nodes in the plane, and hence has an angular resolution π/m. The quadratic scaling of the radial nodes implies that their density decreases from the inner to the outer classical turning point. Therefore, the wave function can detect the smallest spatial features on a circle given by the inner classical turning point, R min . The (angular) resolution corresponds to the distance separating two adjacent nodes on this circle, w = R min sin (π/m). Adapting w to a square lattice gives a critical lattice spacing equal to w/ √ 2. Dropping terms of order ν −2 and assuming the angular resolution is smaller than any radial wave function feature we arrive at a critical lattice spacing The black curves in Fig. 5 are lines through the points (d β,m , E β,m ), where the connection between β and l is made with the approximate relation l = 2(β − 1) − m, as discussed in Sec. IV A. These fit well with the qualitative transitions seen in the spectrum for the lower bands where the approximations are more accurate. The transition can therefore be interpreted as the minimal spacing of scatterers which still can be resolved by the wave function. This spacing d c should not be confused with d D from Table I, which is the maximal spacing for breaking the degeneracy of all levels in the manifold ν. To keep the number of shifted states constant at b D = N D , considerations for the remainder of this section will refer to d ≤ d D .
B. Between the homogeneous and the few-scatterer limit: Chaotic spectra Random matrix theory (RMT) is an appropriate framework to analyze chaotic spectra. Although the Rydberg Composite's spectra are in principle chaotic, the application of RMT to the present problem is hindered by the fact that for d > d c the system obeys several symmetry constraints when F = 1. Moreover, towards the "trilobite-limit" of only a few scatterers (F ≈ 0 and or d d c ) the spectrum becomes regular. Both properties affect strongly the mean density of states. Hence, standard tools [34] from RMT to describe properties of a classically chaotic system are cumbersome to implement as they require knowledge of the mean DoS for unfolding. The unfolded DoS has uniform mean density [35,36] and can be used to extract the eigenvalue correlations in the spectrum.

The adjacent gap ratio (AGR)
To avoid unfolding, we resort to the so called adjacent gap ratio (AGR) [35,37] AGR = min(s n , s n−1 ) max(s n , s n−1 ) with s n = E n − E n−1 and the average taken over the whole spectrum. When F = 1 we also average over many lattice realizations. Since the AGR only depends on local fluctuations it does not require unfolding [35]. AGR can also deal with a "mixed" chaotic and regular spectrum, i.e., it can differentiate Poisson statistics, marking uncorrelated energies typically from preserved subspaces due to symmetries, from Gaussian Orthogonal Ensemble (GOE) statistics which occur for chaotic dynamics [34] without additional symmetries, when level repulsion is present [36]. As RMT references for regular and chaotic Rydberg Composite dynamics we obtain, for a matrix size corresponding to the ν = 30 case, AGR P = 0.386 and AGR GOE = 0.530 for Poisson statistics and GOE statistics, respectively. For further details see Appendix E.
2. The evolution of AGR with the fill factor for different fixed lattice spacings One can see in Fig. 9a that towards small fill factors, but compatible with N > N D where the spectrum looks chaotic, the AGR function g d (F ) indeed approaches g 0 ≡ AGR GOE for all d shown, see Appendix E. However, g d (F ) breaks off g 0 for increasing F , to reach eventually the value g d (1) = 0 due to geometry induced degeneracy. For larger d the break-off occurs at larger F . For large d the AGR function is box-like in shape with a sudden transition to g d (1) = 0. We note here in passing that to a good approximation the family of AGR functions g d (F ) shown follow the form an interesting relation revealing a self-similar property, whose deeper analysis is beyond the scope of this work.

The evolution of AGR for a filled lattice with decreasing lattice spacing
For a filled lattice F = 1, Fig. 9b reveals that neither GOE nor Poisson values match with the statistics observed for any lattice spacing d. For small d towards the homogeneous limit, the AGR approaches zero again due to the geometrically induced degeneracies. However, even for d > d c the AGR settles to a value different from the GOE one due to the inherent symmetries of our system. We have simulated a chaotic system obeying the inherent "crystal" symmetries by a block diagonal GOE matrix with each block representing one irreducible representation of the C 4v point group (see Appendix E). This synthetically obtained AGR (black line in Fig. 9b) agrees well with the AGR of the true spectrum for d > d c .
We may conclude that the spectral fluctuations in the DoS are indicative of a chaotic system with symmetries separating states into none interacting blocks in the F = 1 limit, while for F = 1 these symmetries gradually break until the spectrum is purely chaotic.

A. Isolation of Rydberg manifolds in the presence of scatterers
Rydberg Composites live in the Hilbert space of a single Rydberg manifold ν. This implies that the excitation must be high enough such that interactions with adjacent manifolds is negligible. We have shown that the normalized energy levels of Rydberg Composites in the dense lattice limit scale as ν −5 , ν −11/2 , and ν −6 in 1D, 2D, and 3D, respectively. These energies must be compared, as a function of ν, with the overall spacing between Rydberg manifolds in order to ascertain the isolation of the Rydberg Composite's manifold. The Hellman-Feynman theorem guarantees that the coupling between energy levels increases inversely to their energetic separation, and so we must confirm that the Composite spectrum does not overlap, or even approach, an adjacent Rydberg manifold. The spacing between Rydberg manifolds decreases as ν −3 , and hence the scaling of the un-normalized Rydberg Composite spectra must fall faster than this value. In 1D, the un-normalized spectrum is: since d ∝ ν. In 2D, we have (for the strongest scaling, ν −11/2 ): again, using d ∝ ν. Finally, in 3D, d ∝ ν 4/3 , and so In all three cases the Rydberg Composite's energies decrease faster than the splitting between manifolds as a function of ν, and hence at sufficiently high ν only the states of a single manifold contribute and the Rydberg Composite exists as described.

B. Experimental choice of ν
The requirement for isolation of the Rydberg Composite's manifold ν stretches certain experimental possibilities and therefore, it might be desirable to choose a scatterer species with a smaller scattering length than Rb or Cs, the current standards. For example, sodium (a s (0) ∼ −5) and lithium (a s (0) ∼ −7) have scattering lengths only a third to a half that of Rb [38].
A high ν is also important in order to reach an experimental regime where our additional approximations -constant scattering length and vanishing inter-atomic potentials -are realistic. The model Hamiltonian of Eq. 3a neglects scattering contributions from higher partial waves and the polarization interaction between the Rydberg core and the scatterers. All of those approximations become more accurate at higher ν. The most serious obstacle is probably the current experimental capability in creating small lattice spacings. If the experiment was performed in an optical lattice this would require quite a large principal quantum number, as the current minimum lattice spacing is around d = λ/6 ∼ 2500a 0 [39]. However, it may be possible to observe experimentally the onset of the chaotic behavior of the Rydberg Composite spectrum as discussed in section V B when shrinking the lattice spacing as far as possible. This chaotic behaviour continues to larger lattices spacings than d D . One can see in Fig. 10 that at the ratio d/ν = 10 the chaotic AGR value is still present; at the minimum optical lattice spacing mentioned above this would require ν ∼ 250 to observe. This is feasible given that Rydberg states with ν ∼ 300 − 500 have been produced [40].

VII. CONCLUSIONS AND FUTURE WORK
In this article we have introduced Rydberg Composites built by coupling a Rydberg atom to a dense distribution of many neutral atoms immersed within the Rydberg wave function. Rydberg Composites provide a systematic interpolation from the trilobite and polyatomic few-body regime to a dense environment with a homogeneous density of scatterers as the asymptotic limit. Rydberg Composites, particularly the 2D monolayer case emphasized here, are a new form of matter intersecting few-body atomic Rydberg physics, quantum dynamics involving optical lattices, and few-body quasiparticle examples from solid state physics.
One can imagine many immediate possibilities to extend this concept. These include more refined geometries of the embedding environment to tune the Rydberg Com-posite spectrum, a goal which is traditionally reached by applying external electric or magnetic fields. Also, localization and decoherence studies are feasible by removing the frozen gas restriction and either shaking the lattice explicitly or allowing it to move randomly at some finite temperature.
Appendix A: Matrix representation in the "trilobite basis" In the context of polyatomic Rydberg molecules it has previously proven useful to perform a change of basis from the manifold of Rydberg states |νlm to the basis of trilobite dimer states |i [41]. Each element of this nonorthogonal basis is a trilobite wave function extending from the Rydberg core to the ith scatterer, This basis is ideal when ν 2 M 1, as it greatly reduces the numerical challenges associated with the large ν 2 -dimensioned Rydberg basis. It furthermore provides significant qualitative insight into the structure and possible geometries of polymers since the eigenstates give directly the contribution of each scatterer within the configuration to that energy configuration. From this one can define alternative localization measures utilizing the information immediately available from these eigenvectors, or classify the states within this basis using the known symmetries of the scatterer configurations, as we do in Appendix C.
In this appendix we extend this method to the Rydberg Composite. Since this system typically has M ν 2 , the trilobite basis is no longer numerically beneficial. It can still provide useful qualitative insight, and in the infinite scatterer limit it leads to an alternative method, solving an integral equation, to obtain the spectrum.
Within our stated approximations, the representation of H in the trilobite basis is the M × M matrix H ii = Υ( R i , R i ). One numerical advantage of this approach is that the matrix element H ii can be expressed using only u n0 (R) and u n0 (R), eliminating the need to evaluate many high-l wave functions when ν 1 [19]. If M > B D the diagonalization of H ii in this representation yields M − B D vanishing eigenvalues in addition to the B D shifted eigenvalues, and one numerical disadvantage lies in distinguishing these from real, but small, eigenvalues. As M → ∞, the dimension of H ii becomes infinite, and hence the eigenvalue equation becomes an integral equation, where V is the scatterer volume in dimension D. Since Υ( R, r) is separable, Eq. A2 has solutions whenẼ is obtained from the determinantal equation This equation can be solved via a numerical root finder.

Appendix B: Parabolic coordinates
The hydrogen atom separates in many coordinate systems, and one should choose a coordinate system that is, if possible, adapted to the geometry of the scatterer distribution. For example, the eigenstate for a single scatterer is nearly proportional to a Rydberg basis wave function in ellipsoidal coordinates [42], and the scatterer operator in the dense lattice limit for 1D and 3D-Composites commutes with the Hamiltonian in spherical coordinates. In the 2D case, the spherical wave functions are clearly not well adapted to the scatterer distribution. Although cylindrical coordinates are well-suited to the 2D-Composite scatterer distribution, the Coulomb potential does not separate in these coordinates. It does, however, separate in parabolic coordinates, These treat parabolas ξ and η on either side of the z = 0 plane democratically, and therefore could be more closely adapted to the 2D-Composite. As the following shows, the Hamiltonian in this coordinate system still must be numerically diagonalized, although it does have closedform analytic matrix elements using the analytic forms for the hydrogen wave functions in parabolic coordinates. We therefore present this calculation not for its direct usefulness to the problem at hand, but to define a potentially useful yet uncommonly employed starting point that could benefit future calculations of the properties of Rydberg Composites. The hydrogen wave function in these coordinates is where n = n 1 +n 2 +|m|+1 and β = 1 n . This is normalized with respect to the volume element, (ξ + η)/4dξdηdφ. The matrix elements of the scatterer potential are Integration over ϕ again leads requires m = m , while integration of ξ sets ξ = η. The resulting expression involves only an integral over η, The allowed quantum numbers are also restricted such that n 1 + n 2 = n 1 + n 2 since m is conserved and n is the same within a single manifold expansion. The general integral of this form has a closed form solution, and β = 2m + 1, α = 2/(n 1 + n 2 + m + 1), and λ = α/2. Thus we obtain (2m+1)×(2m+1) block diagonal matrices since n 1 , n 2 , are related to n and m. Diagonalization of these matrices then yields the spectrum computed in the text. Since the basis was not restricted to reject wave functions with no amplitude in the plane from the begin-ning, as we did in spherical coordinates, ν(ν − 1)/2 of these eigenvalues vanish.

Appendix C: Symmetry adapted orbitals
In this appendix we briefly review the use of the projection operator method which, in conjunction with the trilobite basis representation developed in Appendix A, can be used to obtain the Rydberg Composite spectrum when the scatterer configuration is a member of a molecular point group. The particular utility of this approach is that it leads to a classification of the resulting degeneracies and level crossings in the spectrum in the finite lattice-size regime. The description here follows Ref. [19] and is valid only for the s-wave (contact potential) interactions used here; generalization to p-wave interactions requires additional complications [19]. We obtained the symmetry-adapted eigenstates via the following process: • Identify the relevant molecular point group. For example, the planar square lattice satisfies the C 4v molecular point group.
• Construct the labelled basis of trilobite-functions, v, where v k = Υ( R k , r).
• Every symmetry operator in the point group corresponds to a rotation/reflection matrix, denoted r i . This operator acts on the position vector R k of each trilobite function in v, changing it to a different position vector, i.e. r i R k = R j .
• With this information, define an operator R i which acts not on the position vectors R k but rather on the basis vector v. Its elements are (R i ) jj = δ jk δ j k , where k and k are related by r i R k = R k .
• Eq. 20 of Ref. [19], in conjunction with the point group's character table, yields the projection operatorsP j . These are (when properly rankreduced) M j × M -dimensioned matrices, where M j = Tr(P j ). These traces satisfy j M j = M , and thus describe how the total number of eigenstates are partitioned into each irreducible representation.
• These projection operators are then used to partition the Hamiltonian H ii into block-diagonal form, where each block H j kk is the reduced Hamiltonian for the jth irreducible representation. This is done via the transfomration Finally, each H j kk is diagonalized. The eigenstates of a given j exhibit avoided crossings when a parameter, such as d, changes, while the eigenstates corresponding to different irreducible representations (different j values) exhibit real crossings. To make this concrete, we see for the C 4 symmetry of the plane that exactly half (neglecting "round-off" errors due to the mismatch between lattice points in the square and the circular Rydberg orbit) of the eigenstates are in the 2D E irreducible representation, while the remaining 50% of the eigenstates are approximately evenly split among the remaining irreducible representations, A 1 , A 2 , B 1 , and B 2 .
Appendix D: Further details on the smooth DoS This appendix describes additional technical details regarding the density of states calculated in Sec. IV. We begin with the full expression for the density of states used to make Fig. 8a and b, In this formula, F (x; σ, x i ) is a function to convolve the discrete line spectra with a finite width distribution. In Fig. 8 a Gaussian function was used. This has unit normalization, width σ, and peaks about its mean, x i . Other functions, e.g. box functions, could be chosen as well. As discussed in the text, there are different scaling laws for the width σ and energy levels x i for the different regions -"band" and "overlap" -of the density of states. Specifically, in the "band" region this is handled by setting f = 11/2, g = 1, and b = 1. We found that σ = 0.1 sufficed to achieve the smooth resolution of Fig. 8a. The integrated density of states is In the "overlap" region we set f = 6 and g = 1/2; this rescaling of the widths is necessary since the overlap states are denser. The integrated DOS in this case is One inelegant technical detail stems from the fact that the band and overlap regions span very different energy ranges due to the difference between the ν 6 and ν 11/2 scale factors. As a result, we must apply a global compression of the "overlap" energies by multiplying by a somewhat arbitrary factor, b, and afterwards normalize the amplitude of the overall expression with a factor b −2 . We find that b = 0.1 sets, for this range of ν, the two scaled DOS to lie between the same ordinate and abscissa limits. The fully universal scaling of the whole DOS is accomplished by making b, f , and g functions of E: δN δ E = 1 (D5) For each of these fit functions we have found that a tanh function is sufficient to interpolate between band and overlap regions. spanning the range from v 1 to v 2 with a width w and center x 0 provides a smooth interpolating function to transition between these two regions once these parameters are fit to the data.

Appendix E: AGR
The random matrix AGR values were all calculated by diagonalizing 2000 realizations of real symmetric matrices whose matrix elements were randomly sampled from a normal distribution. We observed the AGR values have a weak dependence on matrix size and hence calculated the values used in this paper on matrices of size 465, corresponding to the ν = 30 case for which the majority of our numerical data was calculated.
The AGR value for the GOE case was found to be AGR GOE = 0.5304 ± 0.0003 using dense random matrices. The Poisson value was found to be AGR P = 0.3864 ± 0.0003 using a random matrix with elements only down the diagonal.
In the case of F ≈ 1 the symmetries of the system need to be taken into account. There are two limiting cases, the homogeneous and large d case. The AGR value in the homogeneous case is trivially zero since each case is doubly degenerate due to the ±m symmetry. The system in the large d case belongs to the C 4v symmetry group. The C 4v character table has 5 different types of irreducible representation in it and only states in the same irreducible representation can interact with one another. One of the irreducible representations is two-dimensional, meaning that each state belonging to it is doubly degenerate.
The AGR value for the large d case is calculated using matrices constructed of six GOE matrices in a block diagonal format. The size of each block is chosen to match the number of states in each symmetry irreducible representation found in Appendix C. Four of the blocks are of size 56 (to account for the one dimensional irreducible representations) while the last two are identical (to account for the two-dimensional irreducible representation) and of size 120 each. From this we obtain a value of AGR d large = 0.1894 ± 0.0002.