Emergent percolation length and localization in random elastic networks

We study, theoretically and numerically, a minimal model for phonons in a disordered system. For sufficient disorder, the vibrational modes of this classical system can become Anderson localized, yet this problem has received significantly less attention than its electronic counterpart. We find rich behavior in the localization properties of the phonons as a function of the density, frequency and the spatial dimension. We use a percolation analysis to argue for a Debye spectrum at low frequencies for dimensions higher than one, and for a localization/delocalization transition (at a critical frequency) above two dimensions. We show that in contrast to the behavior in electronic systems, the transition exists for arbitrarily large disorder, albeit with an exponentially small critical frequency. The structure of the modes reflects a divergent percolation length that arises from the disorder in the springs without being explicitly present in the definition of our model. Within the percolation approach we calculate the speed-of-sound of the delocalized modes (phonons), which we corroborate with numerics. We find the critical frequency of the localization transition at a given density, and find good agreement of these predictions with numerical results using a recursive Green function method adapted for this problem. The connection of our results to recent experiments on amorphous solids are discussed.


I. INTRODUCTION
Despite decades of experimental and theoretical efforts, the mechanical properties of disordered solids remain poorly understood. Long-standing challenges such as explaining the low-temperature properties of glasses and ultrasound propagation in granular media require a robust understanding of the effect of disorder on vibrational modes and their localization properties. Interest in phonon localization in disordered solids has been rekindled by recent experiments on colloidal glasses that mimic several vibrational properties of molecular glasses on larger length and time scales [1][2][3]. A significant advantage of colloidal glasses is that their vibrational modes can be experimentally reconstructed, and they are found to exhibit fascinating localization properties. Furthermore, the vibrational modes have been shown to be relevant to understanding the irreversible rearrangements of the system [4][5][6][7][8].
Since Anderson's seminal work [10], the localization of electronic wave-functions in disordered metals has been the subject of several successful theoretical investigations [11]. The study of the vibrational modes of a disordered network of masses and springs (see Fig. 1 ), which is analogous to the electronic problem, received significantly less attention. Formally, the spring problem is different from the electrons problem due to mass conservation, leading to the presence of a zero energy mode -see Eq.
Many important questions remain open regarding the nature of the localization transition of the vibrational modes in disordered systems, its dependence on the dimensionality and the type of the disorder (non-uniform masses versus non-uniform springs, for example). Most studies rely on numerical solutions or study a particular asymptotic limit, and analytic solutions for tunable disorder and system dimensionality are scarce. In this work, we present a model that is analytically tractable, leading to insights into the way the localization transition depends on the disorder magnitude and dimensionality, and elucidating the emergence of a disorder-dependent length scale which determines various elastic properties of the system. Continuum elastic theories of disordered solids typically break down below a mesoscopic length scale * (controlled by the amount of disorder) intermediate between particle diameter and system size [45][46][47]. In jammed packings of grains and random networks, * exhibits a power law divergence whenever the network connectivity with a network of identical masses placed randomly, connected by springs with broadly distributed spring constants -their value exponential in the Euclidean distance between each pair of masses. Only springs with a strength exceeding a threshold Kopt, defined in the main text, are drawn as thin lines. (b) A percolating network of connected strong springs is obtained, yielding a meshwork of pathways, with characteristic length scale l * [9]. (c) The disordered network is coarse-grained to an effective ordered elastic network, with masses m * , spring constants Kopt and length scale of order l * . For system size L l * and dimensions d > 1, the lowlying modes follow continuum elasticity, leading to a Debye spectrum.
or packing density are lowered towards the critical values necessary to hold these systems rigid. At low densities (weak connectivity) the effect of disorder is stronger. Plane-wave vibrational modes, which satisfy the Debye density of states, exist only below a critical frequency ω * ∼ 1/ * that goes to zero at the critical point.
We present a minimal model for a random elastic network that is able to capture these properties. The model is presented in detail in section II. In section III, we argue for Debye spectrum at low frequencies, for dimensions higher than one. Sections IV and V are the core of this work, where we use percolation theory, as shown schematically in Fig. 1, to analytically study the nature of the eigenmodes in the case of strong disorder. We find that above two dimensions there is always a localization transition at finite frequency, for arbitrarily large disorder. We find an approximate analytical expression for the critical frequency, which we verify numerically using a recursive Green function and finite-size scaling technique, presented in section VI. For the delocalized modes, we calculate, analytically and numerically, the speedof-sound, which is directly related to the low-frequency spectrum. The mesoscopic length scale * emerges naturally from the percolation picture, and plays a significant role in the calculation of the speed-of-sound and the phase diagram. In section VII, we discuss the implications of these findings for experiments on colloidal glasses and other amorphous solids.

II. THE MODEL
Our model of the vibrational properties of disordered solids is schematically illustrated in Fig. 1a, which shows point particles (with equal masses, m = 1) randomly and uniformly distributed in a d-dimensional space, with average nearest-neighbor distance r nn . Between every pair of points {i, j} there is a spring with spring constant K ij . We choose: where r ij is the distance between particles i and j, U is a characteristic energy scale and ξ is the range of the exponential interaction. The vibrational frequencies ω i are related to the eigenvalues λ i of the N × N matrix K, as ω 2 i = −λ i . This model describes scalar elasticity, i.e., it neglects the vectorial nature of the forces [24,26]. The randomly chosen points represent the equilibrium locations for the masses. The dimensionless ratio ε ≡ ξ/r nn controls the strength of the disorder -the low density limit (ε 1) corresponds to strong disorder. The reduction to scalar elasticity is not a controlled approximation, but rather a commonly used simplification [48,49]. Nonetheless, previous works suggest that the essential properties of the system are preserved; for example the boson peak and its scaling properties can be studied within scalar elasticity [50], while the flat diffusivity leading to a plateau in the thermal conductivity was first predicted within a scalar elastic model [51].
The sum of every row and column of K vanishes, due to the negative diagonal elements. This 'sum-rule' is a result of momentum conservation, which has important consequences for the localization properties and the resulting vibrational spectrum [52,53]. For example, this constraint makes the matrix negative-semi-definite, which guarantees the system vibrates around a stable equilibrium and leads to a purely delocalized mode with ω = 0, corresponding to a uniform translation of all masses. Without the sum rule, this model is the Lifshitz model for energies in disordered semiconductors [54]. It is precisely the difference between the electronic and phonon problem that enables us to use the percolation formalism here. As we will discuss, the sum-rule makes the problem mathematically analogous to that of electrons incoherently hopping between sites, for which the percolation approach has been successfully applied [55].
Anderson localization is a wave phenomenon, and as such is expected to occur for the classical phonons studied here, for sufficiently strong disorder. Yet it is not obvious a priori whether the low-frequency vibrational modes of this system should be localized or delocalized; It was previously shown that the states with ω > ω * are localized, with ω * → 0 at low density [24]. However, since there is always a delocalized mode at ω = 0, from con-tinuity one might expect that in a large enough system, modes of sufficiently low frequency would also be delocalized. We shall show that this is indeed the case above two dimensions. Another, perhaps simpler, question regards the structure of the density-of-states (DOS) of the eigenvalue distribution at low frequencies, related to the diffusive properties of the associated random-walk problem [56]. Will it be the same as that of an ordered system, manifesting a Debye spectrum, P (ω) ∼ ω d−1 c d (with c the speed-of-sound), or can the disorder change it qualitatively? One may think that a Debye spectrum and delocalized modes should go hand-in-hand, but we shall shortly give a counter-example in a one-dimensional system. In dimensions higher than one we show that one always obtains a Debye spectrum at sufficiently low frequencies.
In one dimension and at low densities, the model we study reduces to a chain of equal masses connected by random springs. For low densities we can neglect matrix elements beyond nearest-neighbors, and obtain a powerlaw distribution of the remaining ones [14,15]. For sufficiently small disorder (large enough ε), one obtains a Van-Hove singularity of P (ω) = const at small eigenvalues, which is Debye-like, while for densities lower than the critical density ε = 1 a sharper power law P (ω) ∝ 1/ω α occurs, with 0 < α < 1. In the one-dimensional case, all the eigenmodes are localized for any density, with the localization length diverging as ω → 0. The fact that the spectrum can still be Debye-like for ε > 1 shows that one can have localized states and Debye simultaneously. Note, however, that at high densities the nearestneighbor model and the one discussed here do not match, since in our case far neighbors will also be important. We now proceed to analyze the model in higher dimensions, focusing mainly on low densities.

III. DENSITY-OF-STATES AT LOW FREQUENCIES
We argue that the density-of-states corresponds to a Debye spectrum at low enough frequencies and for d > 1, for arbitrarily small ε. To show this, it will be useful to employ two mappings between the vibrational problem studied here and 1) a classical random walker in a disordered landscape [41,57] and 2) electrical networks [55].
In the random-walker problem, K describes the Markov process, and the probability of being on site j satisfiesṗ j = K ji p i ; the sum-rule now enforces probability conservation. For this diffusion problem, the return probability ψ i (t) is defined as the probability to be at site i at time t given that the walker was at site i at time t = 0. It can be shown that the return probability scales as [24] is the average over all sites i, and P (λ) is the density of states of the matrix K. This formula implies that a higher DOS at small frequencies corresponds to a larger return probability at long times, which suggests slower diffusion. In general, a Debye spectrum (in any dimension) corresponds to the case of normal diffusion and the return probability scales as 1/(Dt) d/2 . Thus, to establish the existence of a Debye spectrum at low frequencies, it suffices to prove that at asymptotically long times a particle will undergo normal diffusion.
To perform this step, we employ the second mapping, between random walks and electrical networks. Consider a random resistor network whose resistance between two sites is proportional to the matrix elements of K ij , i.e., the hopping rate. Then, the Einstein relation tells us that the diffusion coefficient of the network (if it is non-zero) will be proportional to the conductivity. Therefore, the question of establishing normal diffusion is equivalent to proving that a finite conductivity exists in the network. This result was proved using percolation methods [55], in an identical system where the matrix elements (resistors) between two sites depend on both their energy difference and distance.
For the case of very high (infinite) temperature, the matrix elements depend only on the distance, exponentially, as in the case we discuss here. Thus, for dimensions higher than one, the result of Ref. [55] shows that there is conductivity, and therefore diffusion, at asymptotic times. Moreover, this allows us to derive the diffusion coefficient, since percolation theory gives the dependence of the conductivity on ε: σ ∝ D ∝ e −P d /ε [55].
We now explain how this low-frequency Debyespectrum relates to the results of Ref. [24], where a strong (non-Debye) peak was shown to occur in the DOS at very small frequencies. There, it was shown that for low densities (ε 1) the DOS can be approximated by where V d is the volume of the d-dimensional sphere. This result was derived by approximately calculating all moments of the eigenvalue distribution to an accuracy of O(ε m ) with m > 1. If the correct DOS is a Debye spectrum at low frequency and Eq. (2) for ω > ω , then for sufficiently small ω the corrections to the moments will be sufficiently small to be consistent with the previous result. We expect that ω ≈ ω * , where ω * , the delocalization transition frequency described in section V, is sufficiently small. Although we have shown that above one dimension we expect Debye behavior to set in at sufficiently low frequencies, the existence of Debye spectrum does not prove that the eigenmodes are delocalized; indeed, the onedimensional example in section II, at low disorder, has all eigenmodes localized while the low frequency spectrum is of the Debye form. What happens in higher dimensions? In Appendix A, we show that for the case of small disorder ε 1, the eigenmode behavior is reminiscent of ordered systems: plane waves are approximate solutions to the eigenvalue problem [26] for | k|r nn 1 and ε 1, where k is the wavevector. In Appendix A we also calculate the speed-of-sound of these delocalized modes, a result which agrees well with exact numerical diagonalization, as we show in Fig. 2. Yet it is not clear whether the delocalized nature of the eigenmodes at low frequencies should carry over to the highly disordered case of ε 1 -clearly the eigenmodes cannot be approximately plane waves as they are in the low disordered case. Nevertheless we now use a coarse-graining argument to claim that while locally the modes are very different from plane waves, they resemble plane waves when averaged spatially. We use percolation theory to calculate the length scale over which one has to average to achieve this coarse-grained result.

IV. STRUCTURE OF THE EIGENMODES: COARSE GRAINING AND PERCOLATION APPROACH
To determine the low-energy, large-scale modes of the system, we study the response to a static force applied to the outside of a large but finite system; that is, we consider the compressibility. We follow the method for the equivalent problem of random resistor networks [9], in which the conductivity corresponds to our model's inverse compressibility. The compressibility of the full system can be well described by considering only a subnetwork of springs, where we keep only K ij ≥ K and set all the other K ij to zero; see Fig. 1. The full network, described by K = 0, includes many weak springs, which do not contribute to the rigidity of the system. If K is large, the subnetwork will not connect the edges of the system, so the compressibility will be infinite. The subnetwork just connects the edges of the system when K = K c = e −rc , where r c is the critical radius of the d-dimensional random-site percolation problem [9,55], and η c is the standard continuum percolation threshold [58,59], giving P 2 = 1.20 and P 3 = 0.87, and we have taken ξ = U = 1. This subnetwork is still not a good estimate for the full network's compressibility, because the external pressure is concentrated into a large force at a small number of critical links. Reducing K to a value smaller than K c gives a large number of parallel rigidity paths, each of which has a critical spring with spring constant approximately equal to K c . This results in a mesh qualitatively represented in Fig. 1b. The full network with K = 0 will have a compressibility close to that of this critical subnetwork with K = αK c ≡ K opt , for α a constant of order unity, independent of ε [9, 60]. Reducing K further does not significantly change the compressibility, because the weak springs added are 'shunted' by the critical subnetwork.
The critical subnetwork has a natural length scale giving the size of its "holes" (see Fig. 1b): the percolation length l * = Aε −1 |(r opt − r c )/r c | −ν , for r opt = − log K opt close to r c [9], where A is a constant of order unity and the critical exponent ν is 4/3 in 2D [61] and 0.875 in 3D (derived from [62]). Since α is O(1), |r opt − r c | = | log α| is small, justifying the use of the power-law form: l * ∼ ε −(ν+1) , which diverges at low density.
All subsystems much larger than l * have properties (e.g., compressibility) similar to the infinite system, so we can coarse grain the system into boxes with a size of order l * , as illustrated in Fig. 1c. For determining the system compressibility, each box can be approximated as a spring with a spring constant close to K opt , corresponding to the critical spring in the network connecting the edges of the box. Due to the exponential spread in K ij , all other springs are effectively stiff (K K opt ) or broken (K K opt ) [9,55]. By coarse-graining over regions of order l * , the resulting distribution of critical springs will be narrow, and we have essentially replaced the disordered spring network by an approximately ordered one.
This procedure is only defined in dimensions greater than one, where percolation occurs. For phonons with wavelength large compared to l * , the disorder is averaged away, and the problem maps to that of a weakly disordered system; this is the mechanism by which delocalized modes may come about. The lowest non-trivial eigenmodes resemble plane-waves when coarse grained at a scale larger than l * , as shown in Fig. 4 of Appendix B. At smaller lengthscales, however, they are highly affected by the local environment, and appear disordered. At higher frequencies (but still below the localization transition) modes can be hybrids of plane-wave-like mode [52]. In the following we use the coarse-graining framework to calculate the speed-of-sound of the plane-wave-like modes and the position of the localization transition.

V. SPEED-OF-SOUND AND LOCALIZATION TRANSITION
Using the coarse-graining construction, we find an ordered system on length scales l * , and thus we expect to find large-scale delocalized plane-wave modes. We can readily find the speed-of-sound of these modes at low densities, where percolation applies: c = l * K opt /m * , where m * is the effective mass of each coarse-grained oscillator. Since the volume of each coarse-grained box is comparable to l * d , we expect the total mass of each block to scale as m * ∼ ε d l * Dp ∼ ε d−(ν+1)Dp , where the fractal dimension of the percolating cluster D p is smaller than the system dimensionality d, and is approximately 1.9 in 2D and 2.5 in 3D [63].. Combining these results for l * , m * and K opt gives: The framework presented here also also allows us to find where the localization/delocalization transition should occur; The coarse-grained picture cannot describe features smaller than l * , so we expect the delocalized plane-wave modes to cease to exist as their wavelength approaches l * . We can estimate the delocalization transition frequency Taking the Debye spectrum below ω * , P (ω) ∼ ω d−1 c d , and using the value of c from Eq. 3, we can estimate the number of delocalized modes, and find that it scales as ε νd . As one goes to higher disorder (lower ε) the fraction of delocalized modes decreases, becoming negligible at large disorder.

VI. NUMERICAL ANALYSIS OF THE EIGENMODES
We perform numerical studies to test Eqs. 3 and 4. The prediction of Eq. 3 is corroborated numerically in 2d in Fig. 2, comparing the exact diagonalization of large matrices. For each ε ≤ 0.2, we determined c by averaging results from 2000 realizations with N = 20000 sites with free boundary conditions. For each ε > 0.2 we used a single realization with N = 10000. For each matrix, the low-lying eigenvectors were determined. We discarded the uniform mode and Lifshitz modes [54], where a small number of sites are isolated from their neighbors. The smallest remaining eigenvalue −ω 2 0 gave c = Lω 0 /π, where the system has linear size L = N 1/d /ε, taking ξ = 1. For ε from 0.2 to 0.044, the standard deviation in values of c ranged from 1% to 30% of the mean, respectively. The points shown in Fig. 2 are well determined, with statistical error smaller than the symbol size. Appendix B shows further results that the degeneracies of the low-lying eigenvalues match what is expected from plane waves in an ordered box geometry, see Table S1. Moreover, Fig. 4 in Appendix B shows the average spatial profile of these low frequency modes, which is planewave-like.
Following the same reasoning as the standard one for ordered systems, the DOS which follows from these plane-wave-like modes will be a Debye spectrum, and at low frequencies we find that P (ω) ∼ ω d−1 c d . This is consistent with the results of the alternative derivation for the form of the DOS at low frequencies presented in section III, which relies on mapping this problem to the diffusion problem of a classical random walker in a random landscape. We have used a recursive Green function and finite-size scaling technique to systematically investigate whether states are localized or delocalized in the infinite-size system [64,65], with results shown in Fig. 3. We adapted the technique from the one used in Ref. [ 66], which studied the same system without the diagonal terms of K; see Appendix C for details. The percolation prediction of Eq. 4 indicates that the delocalized modes should persist to arbitrarily small ε. Due to numerical precision limitations, we cannot study ε ≤ 0.05 in three dimensions and thus cannot numerically determine if there is a smaller ε at which all states are localized. Nonetheless, the prediction of Eq. 4 is in good agreement with the numerically determined boundary, as shown in Fig. 3, for both two and three dimensions. The numerical results for 2d also suggest a phase transition, in contrast to the behavior of electronic systems, where it is known that arbitrary small fluctuations will lead to localization in the orthogonal ensemble (which is the appropriate ensemble considering the symmetries of our problem) [67]. Further numerical and analytical study is needed, however, in order to disprove or establish the existence of delocalized eigenmodes in d = 2.

VII. CONNECTION TO AMORPHOUS SOLIDS AND COLLOIDAL GLASSES
Despite its simplicity, our model captures well various features of realistic amorphous solids, such as (a) the transition between Debye and non-Debye behavior in the spectrum (b) the corresponding change in the na- ture of the eigenmodes (delocalized versus localized) and (c) the dependence of the speed-of-sound of the delocalized modes on the connectivity of the system. For example, the low frequencies regime of Debye spectrum has proved very difficult to observe experimentally in colloidal glasses [68]. Our work shows analytically that the extent of this regime can be extremely small for strongly disordered systems but it exists nonetheless. Our parameter ε describes the degree of connectivity of the system, and our predictions for the dependence of ω * , the crossover frequency to a Debye spectrum, qualitatively explains the dependence observed in amorphous solids, where (Z − Z c ) plays the role of ε (with Z the coordination number).
The emergence of l * from our theory is another important feature which has been found in real amorphous solids, but never connected to percolation theory. Both in experiments and within our model, there is a powerlaw divergence of l * as a function of the 'sparsity' of the system (described by (Z − Z c ) or by ε). In both cases ω * ∼ 1/l * . We show that within our model the speedof-sound vanishes as ε → 0, with direct parallel to the "softening" experimentally and numerically observed in jammed solids close to Z c [69][70][71]. Therefore, our model provides a natural framework to predict the mechanical properties of these fragile solids close to the threshold of vanishing rigidity.
It should be emphasized, though, that this is a simplified model, amenable to analytic treatment -it should not be expected to predict quantitatively the transition frequencies in real amorphous solids. Rather, it illustrates that there should be such a transition in the character of the vibrational modes and how it depends on the dimensionality and connectivity of the system. For instance, our analytic framework predicts the dependence of the localization properties on dimensionality, showing that d = 2 is a critical dimension, and that the behavior in 1d is very different from higher dimensions, because there is no percolation in one dimensional systems.
We believe that further insights in the physics of amorphous solids can be gained from studying further this minimal model and compare its predictions for the thermal conductivity and phonon scattering with relevant experimental studies of glasses at low temperatures.

VIII. CONCLUSIONS
We have studied the vibrational spectrum of a disordered system, using analytical arguments relying on percolation theory and a coarse-graining analysis, numerically exact diagonalization, and finite-size scaling. We find that for d > 1, there is a critical frequency ω * vanishing with density, below which a Deybe DOS is observed. This frequency marks a localization-delocalization transition of the vibrational modes. In two dimensions further studies are needed in order to ascertain whether this is a true phase transition. We analytically account for the speed-of-sound of the delocalized modes, both for the low and high density regimes, using a percolation approach. A lengthscale l * emerges, which affects the speed-ofsound and the phase transition. This lengthscale appears in various experiments on amorphous solids, and we discuss the applicability of our work to such systems. In the future, it would be interesting to study heat transport in these systems and to compare results of the minimal model presented here with those of amorphous solids. * These authors contributed equally to this work.

APPENDIX A: SPEED-OF-SOUND AT HIGH DENSITIES
At high densities, Ref. [26] gives a simple argument showing that plane waves solutions can be found at high densities. Indeed, the eigenvalue equation is j K ij u j = λu i , which can be written as: where f (r) = exp(−r/ξ). Using the fact that K ii = − j =i f (r ij ), and guessing an eigenmode u j = e i k· rj , the eigenvalue equation becomes Provided that the sum can be well approximated by an integral, the LHS is independent of i, showing that it is indeed an approximate eigenmode. There are two conditions for this to be valid: 1) The points should be dense enough such that the function f (r) is well sampled (i.e., that we can view the sum as a discretized version of the integral). This is equivalent to the condition ξ r nn , i.e., ε 1. 2) In a similar fashion, the points should be dense enough that the function e i k· r is well sampled. This leads to the condition | k| 1/r nn . If conditions (1) and (2) are met, we find that: If we denote the d-dimensional Fourier transform of f ( r) by F ( k), this leads to [26]: For our exponentially decaying matrix elements, in 2D we find that for long wavelength phonons with 1/k ξ, ω = ck, with: where U = ξ = 1. Thus, these long wavelength phonons are indeed acoustic ones, with a well-defined speed of sound c. This result is verified numerically in Fig. 2.
Repeating the same procedure for 3D, we find:

APPENDIX B: NUMERICAL VERIFICATION OF DELOCALIZED MODES
We use two numerical techniques to verify the existence of delocalized modes in this model: diagonalization of finite systems and a recursive Green function method with finite size scaling. In this section, we describe the former calculations, and in the next section we describe the latter. In both cases, the numerical analysis must truncate Eq. 1, so the elements K ij are set to zero for r ij > r cutof f .
For sufficiently large cutoff values (generally 3-5 times r c , the critical percolation radius, described in section ), the results are independent of the value of the cutoff, as they should be.
For the exact diagonalization, we must discard some low-energy modes. When a single point is further than r cutof f from its nearest neighbor, the resulting eigenspectrum of K contains a mode perfectly localized at that point with frequency zero. Without r cutof f , such modes would still exist, with frequencies exponentially small in the distance to other sites; such modes are called Lifshitz modes [54]. When looking for low-lying modes, we remove from consideration all such Lifshitz modes, defined as those with participation ratios of order unity with vanishingly small frequencies; such modes are also overrepresented in calculations with free boundary conditions, which are the most convenient to use in the degeneracy analysis below. The contribution of these modes to the density-of-states is negligible for a large enough system and for small enough frequencies, since it goes to zero faster than any power-law (note that this relies on the exponential decay of the matrix elements, and for other forms decaying faster with the distance this assertion might be false).
Aside from the Lifshitz singularities, the low frequency modes appear to be delocalized in these finite systems in two and three dimensions. In two dimensions, the lowest two non-trivial frequencies are nearly degenerate, and in three dimensions triply degenerate, as expected from the analysis of eigenmodes in an ordered square lattice geometry with free boundaries. Looking at higher nontrivial eigenfrequencies, their values are consistent with the coarse-grained ordered picture, predicting particular ratios of the low-lying frequencies, as shown in Table 1.
The average spatial profile of these low frequency modes clearly shows plane wave oscillations, as shown in Fig. 4.

APPENDIX C: RECURSIVE GREEN FUNCTION CALCULATIONS
The recursive Green function technique with finitesize scaling was first used to study Anderson localization by MacKinnon and Kramer [64]. It studies quasi-onedimensional systems with one very long axis and d − 1 shorter axes of width w. The system is built recursively by stitching together many (d − 1)-dimensional strips. In the original Anderson model, each site is coupled only to its nearest neighbors, so each new slice has direct couplings only to the closest neighbors; this property is essential for the recursion.
In Ref. [72], we extended this technique to non-lattice problems, as studied in this paper. Without a lattice, the Index ω, 2D ord. ω, 2D num. ω, 3d ord. ω, 3d num.  4. The lowest frequency mode is shown for a single matrix with N = 500, 000, ε = 0.1, in two dimensions. The amplitudes were binned along one of the sides of the square in which the points are randomly chosen, and the graph shows the sum of amplitudes in each bin. While the eigenmode itself is not a simple plane wave, after 'coarse-graining' in this way one obtain a cosine profile, precisely as would be obtained in an ordered system.
quasi-one-dimensional system is stitched together from slices of width w and length l, which contain a random number of sites, chosen from a Poisson distribution with a mean density ε d . In order to have particles in each slice interact only with particles in the nearest neighbor slices, we must impose a cutoff on the interaction between sites.
That is, we set exp(−r) → 0 for r ≥ L c for some cutoff L c . Since we use periodic boundary conditions in the d − 1 dimensions of the slice, we choose the cutoff to be L c = min[l, w/2]. Since there is no lattice, w and l must be chosen appropriately large so that the system is not disconnected.
We begin with a single slice of width w and length l.
We add a second slice and find the portion of the energyresolved Green function G(λ) that connects any site in slice 1 to any site in slice 2. We recursively add more sites, always finding the portion of G connecting sites in slice 1 to sites in slice N , where M i is the number of sites in slice i. The probability for a particle injected with energy λ into slice 1 to make it to any site in slice N is Since all states are localized in 1d systems, P N decays like exp(−2N l/ w ), for localization length w . We extract the effective localization length at energy λ and width w by N must be chosen large enough to get good statistical accuracy for w [65,73]. Statistically, the system behaves as though it consists of N l/ w independent and identically distributed samples chosen from a normal distribution, so the sampling error decreases as N −1/2 . When w has been acquired for a range of widths w, the theory of single-parameter scaling is used to extrapolate to the infinite-size system [64,65,74]. According to single-parameter scaling, for sufficiently large w that irrelevant corrections can be neglected, the dimensionless scaling variable Λ w ≡ w /w obeys a universal scaling law for some unknown universal function F , where φ(ε, λ) characterizes the infinite-size system: it is the localization length for localized systems and the correlation length for delocalized systems. In particular, if Λ w increases with w, then the state is delocalized and if Λ w decreases with w, then the state is localized. This is the approach adopted to determine the phase diagram in Figure 3. Single-parameter scaling has been demonstrated to hold in similar systems without the sum rule [72], but it has not yet been proved in this system. It is possible that the emergence of the percolation length l * would cause a breakdown in the single-parameter scaling, but the data so far do not show any such signs. The details of setting up the recursive Green function calculation have been extensively given elsewhere [75,76]. The basic principle is that, in each step, the exact Green function G N of the wire with N slices and the exact Green function of the new slice g are known. We write G 0 = G N + g. The interaction V that connects them can be treated non-perturbatively using a Dyson equation Matrix elements taken on Eq. 13 give the recursion relations required.
For our case, the sum rule in the definition of the matrix K means that the interaction V contains not only off-diagonal terms connecting sites in the two slices, but also contains diagonal terms, adjusting the on-site energy of each site. This does not pose any difficulty for the method, but requires more matrix operations per calculation, increasing the cost of the calculation.
In each calculation, we choose the values of λ we wish to study, and for each λ store the matrices G N 1N , G N 11 and G N N N . For the N + 1st slice being added to the system, it has a Hamiltonian K and the off-diagonal interaction connecting the existing wire to the new slice is U . Then in the usual problem, without the sum rule, the result is where we used g −1 = λ − K. When we include the sum rule, the interaction contains the off-diagonal terms U , the diagonal terms U N acting on the sites in the Nth slice, and the diagonal terms U N +1 acting on the sites of the new slice. Then we find If the calculations are done all at once, there is no need to store G 1,1 . But it is often convenient to break the calculation of a very long wire into many smaller pieces, and then stitch them together. This stitching requires that G 1,1 be stored. Consider one wire with N 1 slices with recursively calculated Green functions G N1 1,1 , G N1 1,N1 , and G N1 N1,N1 and a second wire with N 2 slices with recursively calculated Green functions g N2 1,1 , g N2 1,N2 , and g N2 N2,N2 . Let N = N 1 + N 2 . Then we can stitch g onto G to find G N N,N = g N2 N2,N2 + G N N1+1,N2 U g N2 1,N2 + (g N2 1,N2 ) † + G N N1+1,N2 U g N2 1,1 ) G N N1+1,N2 = g N2 1,N2 1 1 N1+1 − U N1+1 g N2 It is generally desirable to let N be large enough that the statistical error in w is less than 1%. For values of λ, ε where the states are highly delocalized, w can be very long, requiring an inordinate length N to get reasonable errors. Conveniently, the values of w are reasonable for the localized states and the nearby delocalized points, so we can extract the delocalization boundary ω * (ε), while we cannot quite prove that the strongly delocalized states are actually delocalized. Care must be taken to avoid both overflow and underflow in storing G 1,N , as detailed in Ref. [ 76]. Figure 5 shows some of the calculations used to create Figure 3. The statistical error bars are smaller than the data points for all but the lowest energy (most delocalized) state. For these simulations, l = 50 and N ranged from 2 · 10 5 to 2.6 · 10 7 .  Figure 3. For ε = 0.097 in d = 2, we show Λw(w). The phase boundary between the localized (sloping down) and delocalized (sloping up) energies is clear.