Order to disorder in quasiperiodic composites

From quasicrystalline alloys to twisted bilayer graphene, the study of material properties arising from quasiperiodic structure has driven advances in theory and applied science. Here we introduce a class of two-phase composites, structured by deterministic Moiré patterns, and we find that these composites display exotic behavior in their bulk electrical, magnetic, diffusive, thermal, and optical properties. With a slight change in the twist angle, the microstructure goes from periodic to quasiperiodic, and the transport properties switch from those of ordered to randomly disordered materials. This transition is apparent when we distill the relationship between classical transport coefficients and microgeometry into the spectral properties of an operator analogous to the Hamiltonian in quantum physics. We observe this order to disorder transition in terms of band gaps, field localization, and mobility edges analogous to Anderson transitions — even though there are no wave scattering or interference effects at play here. Moiré patterns and ordered aperiodic geometries have received significant attention since their observation in twisted bilayer graphene and quasicrystalline alloys. Here, the authors theoretically demonstrate that the same patterns can govern localization in quasiperiodic metal-dielectric composites and can be used to engineer the resultant optical and transport properties.

I n the late 1980s it was shown that in a composite patterned after a crystal, such as a dielectric material with a periodic lattice of voids, electromagnetic waves of certain frequencies and directions could be prohibited from propagating within the structure 1,2 . This observation established a powerful analogy relating photonic band gaps to electronic band gaps in metals and other condensed matter. Thus solid state physics and Anderson localization was brought to optics [1][2][3][4] , leading to the development of photonic crystals and theories of controlling the flow of light through structured media. The discovery of quasicrystals [5][6][7] demonstrated that geometries with predictable long range order but no periodicity could play an important role in physics and materials science. This led to the development of photonic quasicrystals [8][9][10][11][12][13][14][15][16][17] , with the conceptual framework again provided by the analogy with quantum transport in solid-state physics.
Motivated by these findings and the highly active field of twisted graphene bilayers 18 , with Moiré patterns tuned by the twist angle to take periodic and aperiodic geometries, here we construct a class of deterministic, two-phase Moiré-structured composite materials in two dimensions. This construction enables us to study in several physical settings how classical transport behaves in the transition from periodicity to aperiodicity. Indeed, rather than a governing wave equation like Schrödinger's equation for quantum transport or the classical wave equation for electromagnetic transport 17,[19][20][21] , problems involving electrical conductivity σ, thermal conductivity κ, complex permittivity ϵ in the quasistatic limit, or diffusivity D can all be formulated in terms of the same divergence form second-order elliptic equation (2) below, and do not involve any wave interference or scattering effects. Bulk behavior is analyzed in terms of the Bergman-Milton (or Stieltjes integral) representation, which holds for the effective parameters σ * , κ * , ϵ * , D * , etc. [22][23][24][25] . It involves a spectral measure μ of a self-adjoint operator G, which plays the role of the quantum physics Hamiltonian and depends only on the mixture geometry. In discrete settings, G is a real-symmetric matrix. The measure μ, local electric field E, displacement D = ϵE and current J = σE are all determined by the eigenvalues and eigenvectors of G. One of our main results is that through this spectral distillation and recent results on computing μ 26 and analyzing its behavior with random matrix theory 27 , we establish a powerful analogy between various classical transport processes in periodic and quasiperiodic composites, and quantum transport with localization and band gaps in solid state physics, as was done for optics in photonic crystals and quasicrystals in the scattering regime. We emphasize, however, that our results apply broadly to transport phenomena in settings described by Eq. (2), with no restriction on the length scales in the systems involved, except for the condition imposed on the microstructural scale by the quasistatic assumption that must be satisfied in the context of complex permittivity.

Results and discussion
We find that as the geometry is tuned from periodic to quasiperiodic, the eigenvalues, eigenmodes, profile of ϵ * , and localization properties of E undergo an order-to-disorder transition analogous to the Anderson transition. Our results are described in the (quasistatic) electromagnetic case, but we keep in mind their broad applicability. Spectral measures for periodic systems have sharp resonances that induce dramatic variability in band and absorption characteristics, and in profiles of ϵ * . Regions of extended eigenstates are separated by "mobility edges" of localized states, and E is localized for certain frequencies and extended for others. As the geometry is tuned to aperiodicity, the behavior of μ and ϵ * resembles that of the 2D random percolation model at its threshold, with a regularly distributed mixture of localized and extended eigenstates giving rise to tenuously connected current paths, pronounced spectral endpoint behavior, and Wigner-Dyson eigenvalue statistics with strong level repulsion 27 .
Our investigation here of quasiperiodic media was motivated not only by the findings for random media in ref. 27 , but by much earlier studies which revealed sensitive, discontinuous dependence of bulk transport on the variations in local properties 28,29 . For example, it was found in one dimension with local conductivity σðxÞ ¼ 3 þ cos x þ cos kx, which is periodic for k rational and quasiperiodic for k irrational, that the effective conductivity σ * (k) is discontinuous in k 28 , with 2D examples in ref. 29 . These studies, in turn, were motivated by the discovery of quasicrystals and findings on the spectrum of Hamiltonians with quasiperiodic potentials [30][31][32] .
The spectral characteristics considered here govern the optical properties of nanostructured bimetallic films 33,34 and depositions of nanosized metal particles on thin dielectric substrates [35][36][37][38] , which change as a function of heterogeneous surface structure composition and geometry. This enables tunability of their optical responses for nano-plasmonic device applications [33][34][35][36][37][38] . The long wavelength quasistatic assumption holds in the visible range 39 , and these systems are described macroscopically by the Stieltjes integral representations for ϵ * or σ * . Resonances in μ explain giant surface-enhanced Raman scattering observed in semicontinuous films 34,40,41 , and induce strong fluctuations in E and the dielectric profile of ϵ * , associated with the excitation of collective electronic surface plasmon modes 39 . We numerically explore these phenomena in 2D impedance networks with quasiperiodic microgeometry and discuss our results using Anderson transition interpretations of random matrix theory.
Effective transport in composites with Moiré-structured microgeometry We begin by introducing a class of 2D two-component composites whose microgeometries are based on Moiré patterns, and are tunable to be periodic or aperiodic.
Constructing Moiré patterns. Consider the square bond lattice joining nearest neighbor points in Z 2 , with standard basis vectors e 1 and e 2 , and the scaled rotation transformation T defined for ðx; yÞ 2 R 2 by T : ðx; yÞ 7 ! ða; bÞ ; The mixture geometry of the two phases is determined by the characteristic function χ 1 , taking the value χ 1 = 1 in material phase 1 and zero otherwise, with χ 2 = 1 − χ 1 . The system microgeometry is constructed from the periodic function ψða; bÞ ¼ cosð2πaÞ cosð2πbÞ and the condition χ 1 (x, y) = 1 for all ðx; yÞ 2 R 2 such that ψ(T(x, y)) ≥ ψ 0 , and is zero otherwise. We focus on the value ψ 0 = 0, which generates in the underlying bond lattice a discretized composite microstructure with a fraction p ≈ 1/2 of type one bonds. We do so to compare features of deterministically tuned quasiperiodic systems to those of the random percolation model near the percolation transition p = p c = 1/2 42,43 .
Primitive translation vectors for ψ are t 1 = (1/2, 1/2) and t 2 = (1/2, −1/2). When r and θ are chosen such that T : ðme 1 þ ne 2 Þ 7 ! ðm 0 t 1 þ n 0 t 2 Þ for integer values of m; n; m 0 and n 0 , then χ 1 has a finite period of, at most, K ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , and has infinite period otherwise. The arrangement of r and θ such that K < ∞ is fractal in nature, as shown in Fig. 1. The arrangement of (r, θ) values associated with finite periods is similar to fractal distributions defined in terms of rational numbers on the real line, such as Thomae's function 29 .
An integral representation for effective transport coefficients.
The effective behavior of macroscopic transport in two-phase composite materials is described by homogenized coefficients including electrical and thermal conductivity, diffusivity, complex permittivity, and magnetic permeability. These can all be defined in terms of the same elliptic partial differential equation 25,43 . For complex permittivity in the quasistatic regime, such as the metaldielectric mixtures in visible light discussed above, the system is described locally by with potential ϕ, electric field E = −∇ϕ, displacement D = ϵE, and local complex permittivity ϵ(x, y) taking frequencydependent values ϵ 1 (ω) or ϵ 2 (ω), where 〈E〉 = E 0 and 〈⋅〉 denotes spatial average. The fields E and D satisfy ∇ × E = 0 and ∇⋅D = 0, with ϵ = ϵ 1 χ 1 + ϵ 2 χ 2 . See refs. 24-26 for a "weak" formulation of this problem that rigorously accounts for the discontinuous, and thus non-differentiable nature of the parameter ϵ(x, y) in Eq. (2). The effective permittivity matrix ϵ * can be defined by 〈D〉 = ϵ * 〈E〉 with 〈E〉 = E 0 , where E 0 = E 0 e k for some standard basis vector e k , k = 1, . . . , d, where d is dimension. Equivalently, it can be defined in terms of system energy using hD Á Ei ¼ ϵ Ã kk E 2 0 , where ϵ Ã kk is the kth diagonal coefficient of the matrix ϵ * , which we denote by ϵ Ã ¼ ϵ Ã kk . Thus, the effective parameter characterizes a homogeneous medium immersed in a uniform field E 0 that behaves macroscopically and energetically as does the inhomogeneous composite medium.
A key feature of Eqs. (3) and (4) is that the material parameters in s and the applied field strength E 0 are separated from the geometric complexity of the system, which is encoded in the properties of the spectral measure μ and its moments μ n ¼ R 1 0 λ n dμðλÞ. For example, μ 0 = 〈χ 1 〉 = p, the volume fraction (or area fraction) of medium 1. All of the effective coefficients of the composite material mentioned above are represented by Stieltjes integrals with the same μ 44 .
While the measure μ can include discrete and/or continuous components 25 , it reduces to a weighted sum of Dirac δ-functions δ(λ − λ j ) for media such as laminates, hierarchical coated cylinder and sphere assemblages, and finite RLC impedance networks [22][23][24][25][26] . Here, we investigate effective transport properties of square two-component impedance networks in 2D of size M with periodic and quasiperiodic microgeometry. In this setting, G = χ 1 Γχ 1 is a real-symmetric matrix of size N = 2M 2 , χ 1 is a diagonal matrix with 1's and 0's along the diagonal corresponding to impedance type, and where ∇ is a finite difference matrix representation of the differential operator ∇ 26 . The measure μ is determined by the eigenvalues λ j and eigenvectors v j of N 1 × N 1 submatrices of Γ with rows and columns corresponding to the diagonal components ½χ 1 jj ¼ 1, with 26,27 . In this case, Eqs. (3) and (4) become finite sums with given explicitly in terms of the λ j and eigenvectors v j of G 26 .
In the next section, we compute the spectral measure μ, hence the local fields and the effective complex permittivity ϵ * for the Moiré-structured class of composite materials described by Eq. (1) above. We interpret the frequency dependent behavior of physical quantities such as the phase and amplitude of ϵ * , and localization and intensity of E and D, in terms of spectral properties of μ and Anderson transition interpretations of random matrix theory.

Analysis
The Moiré system introduced above is parameterized by r > 0 and 0 ≤ θ < 2π, which generates a diverse assortment of periodic ("finite period") and quasiperiodic ("infinte period") microgeometries. To numerically calculate mathematical and physical quantities, we consider finite subsets of these systems as RLC impedance networks. Different types of microgeometries in this class are displayed in Fig. 2a with small enough system sizes to resolve the small-scale geometry while still illustrating the large variety in structure, hinting at the geometric richness of our Moiré composites. The bond color indicates the modulus value of E, i.e., |χ 1 E|, calculated via Eq. (6). Since χ 1 D = ϵ 1 χ 1 E these colors also specify displacement values with a change in scale by |ϵ 1 |. We therefore normalize the computed fields to take values in the unit interval.
It was shown in ref. 26 that expressions known in closed form for the 2D percolation model in the infinite volume limit are well approximated by ensemble averages of systems of size ≈70 and by single systems of size ≈200. The Moiré-structured composites studied here can have coherent structures on large length scales. However, we found for a system size of 199 that fluctuations present for smaller systems have essentially stabilized, with numerical results visually identical for larger system sizes, e.g., ≈250. A systematic study of system size dependence of quantities and finite size effects is interesting and useful but beyond the scope of the current manuscript.
Transport behavior in microstructures along a short trajectory in parameter space. In this section, we investigate a small swath of the large parameter space, for r ¼ ffiffiffiffiffi 10 p =3 and θ ¼ arctanð1=3Þ þ ϕ for 0°≤ ϕ ≤ 2°, starting from a short period system. Figures 2b and 3a display examples of this region of parameter space and show that such a small change in the Moiré twist angle θ gives rise to a dramatic transition in composite microgeometryfrom a short period system with orderly field (or current) paths to quasiperiodic systems with disorderly, meandering paths similar to those exhibited by the random percolation model near p = p c .
When the fields are plotted versus s(ω), 0 ≤ Re sðωÞ ≤ 1 with Im sðωÞ ( 1 a frequency dependent localization/delocalization transition of fields is revealed for small values of ϕ ∈ [0, 2], as shown in Figure 2b for ϕ = 1/8 and 1/2. In contrast, the fields for angles closer to ϕ = 2 are more disordered and resemble those in the random percolation model, and are qualitatively similar to the rightmost panel in Figure 2b for all 0<Re s ≤ 1. We investigate these and other phenomena through mathematical and physical quantities such as the spectral measure μ, correlations of its eigenvalues, localization of its eigenvectors, phase and magnitude of ϵ * , localization and intensity of E, etc. A large variety of physical phenomena exhibited by inhomogeneous materials can be described by two-component RLC impedance networks 40 . Each of the two components is created by combining a resistor R, inductor L, and capacitor C in a way that achieves an impedance characteristic of the material being modeled. For example, a Drude-metal/dielectric composite is modeled by R and L in series, in parallel with C for one component, and C for the other 39 , yielding a plasma frequency ω 2 p ¼ 1=LC and relaxation time τ = L/R. As Kirchhoff's network laws are discrete versions of the curl-free and divergence-free conditions on the fields in Eq. (1), these RLC impedance networks really do resemble the continuum composites they are intended to model 39 .
Indeed, the AC response and polarization effects observed in a variety of conductor-dielectric mixtures at low frequencies are modeled by an R − C network, while metal-dielectric composites exhibiting collective electronic modes at higher, optical frequencies such as (surface) plasmon resonances are modeled by an RL − C network 40,41 . The dependence of s(ω) on frequency ω is model specific. For the R − C and RL − C models, 0 ≤ Re s < 1 for 0 ≤ ω < ∞ and δ ≤ Im s < 0, where |δ| can be chosen as small as desired, with Re s ! 1 and Im s ! 0 as ω increases 40,41 . In order to give a model independent description of the phenomena investigated here, we plot s-dependent quantities using 0 ≤ Re s ≤ 1 and Im s ¼ 0:001 fixed. For the sake of discussion, we describe our results in terms of the optical regime for the Drude model for gold/vacuum composites, which roughly corresponds to the interval Re s ⪅ 0.2. The optical regime for other material combinations corresponds to values of Re s throughout the unit interval 45 .
As the frequency changes and s(ω) sweeps across the complex plane, with s(0) = 0, the spectral measure μ, distribution of its eigenvalues, and localization properties of its eigenvectors, shown in Fig. 3, govern the frequency dependence of the phase and magnitude of ϵ * and the intensity and localization of E and D, shown in Fig. 4, according to the formulas in Eq. (6). Keeping these formulas with Im sðωÞ ( 1 in mind, we call resonant frequencies the values of ω where Re sðωÞ % λ j and the masses m j of μ are largest (shown in red in Fig. 3) and/or there's a large density of eigenvalues λ j with moderate to large values of m j .
For the short period system with ϕ = 0 shown in Figure 3a, the spectral measure μ in Figure 3b is comprised of sharply peaked resonances. As ϕ increases and the composite microgeometry becomes quasiperiodic, the resonant frequencies away from ω = 0 (λ = 0) spread out, change frequency locations, and diminish in strength. As ϕ → 2, the resonances in μ continue to spread out until all but the Drude resonance at ω = 0 diminish, and μ and ϵ * begin to resemble those of the random percolation model for p = p c , shown in the rightmost panels of Fig. 3. Departure from short period system geometry is parameterized by ϕ. Composite microgeometry parameterized by r ¼ ffiffiffiffiffi 10 p =3 and θ ¼ arctanð1=3Þ þ ϕ for 0°≤ ϕ ≤ 2°with system size 199. For small values of ϕ, the fields exhibit a frequency dependent transition from localized (loc) to extended (ext). Identical values of ϕ correspond to identical microgeometries, and the differences in the values of |χ 1 E| are solely due to frequency ω dependent material properties for different values of ω. As ϕ → 2, the local fields become similar for all frequencies away from ω = 0, qualitatively resembling the rightmost panel in (b) (as well as that of the percolation model near the percolation threshold p = p c 27 ). The inverse participation ratio (IPR) providess quantitative description of this localization phenomenon.
Resonances in μ have a physical interpretation in terms of relaxation times in the transient response in the R − C model, or in terms of dielectric resonances in the RL − C model 40,41 . The dielectric resonances observed for the RL − C model with percolative geometry have been argued to provide a natural explanation for the anomalous fluctuations of the local electric field E, which are responsible for giant surface-enhanced Raman scattering observed, for example, in semicontinuous metal films 41 . We show that the resonances in μ shown in Fig. 3 give rise to dramatic fluctuations in the amplitude and phase of ϵ * and the intensity of the fields E and D.
The inverse participation ratio (IPR) characterizes vector localization phenomenon. For an N 1 -dimensional unit vector u it is given by IPRðuÞ ¼ ∑ i u 4 i , where u i is the ith component of the vector u, i = 1, …, N 1 , and satisfies IPR(u) = 1 for a completely localized vector with only one non-zero component and IPR(u) = 1/N 1 for a completely extended vector with all components equal in value 27 . For matrices in the Gaussian orthogonal ensemble (GOE), the eigenvectors are quite extended with a mean asymptotic IPR value of IPR GOE = 3/N 1 27 . Figure 3c displays IPR(v j ) for the eigenvectors v j , j = 1, …, N 1 , of G for various values of ϕ, as a function of the eigenvalues λ j . The red dots in Fig. 3b, c for ϕ = 0 and 1/8 show that resonant frequencies correspond either to very extended eigenvectors or "mobility edges" where the values IPR(v j ) have large variability  for small changes in λ j . As ϕ increases and the microgeometry becomes quasiperiodic, the mobility edges diminish as the values IPR(v j ) become more regularly distributed and qualitatively similar for all 0 < Re sðωÞ ≤ 1 away from the Drude peak at ω = 0, as shown in the two rightmost panels of Fig. 3b. As ϕ → 2, the IPR(v j ) resemble those of the random percolation model at its threshold p = p c = 1/2, as shown in the rightmost panel of Fig. 3c.
The frequencies corresponding to resonances of μ and field delocalization are tunable through the quasiperiodic microgeometry via the scale parameter r and Moiré twist angle θ in (1), which is critical to potential engineering applicationsgiven a desired frequency dependence for the profile of ϵ * and field localization, values of r and θ can be selected accordingly. This is illustrated in Fig. 5 which displays the θ-dependence of the eigenvalue density ρ(λ, θ) in (a), the spectral function μ(λ, θ) in (b), the magnitude and phase of the relative effective permittivity ϵ * /ϵ 2 in (c) and (d) and the IPR in (e), with r ¼ ffiffi ffi 5 p =2 fixed. Short period systems are indicated by dark horizontal streaks due to associated isolated resonances in μ, with localized regions of yellow. Figure 5a shows that some of these resonances in μ(λ, θ) are due to resonances in ρ(λ, θ). However, in Figure 5b, the significance of the measure mass becomes apparent, which can diminish eigenvalue resonances or even create resonances in μ(λ, θ) in regions of low eigenvalue densityalso illustrated in the leftmost panel of Figure 3b by the individual eigenvalue λ j ≈ 0.32 with relatively large spectral mass m j ≳ 0.1. The influence of μ on ϵ * /ϵ 2 is striking with resonances and features in |ϵ * /ϵ 2 | following those in μ, and with an antisymmetry in phase(ϵ * /ϵ 2 ) about Re s % 0:5. The IPR values displayed in Fig. 5e again illustrate that resonances in μ are associated either with extremely extended eigenvectors or mobility edges, with large variability in IPR values for a small change in λ. The symmetry ρ(λ) = ρ(1 − λ) well known for the percolation model 26,41 is evident in Fig. 5a for quasiperiodic geometry, and also has symmetry for θ between π/8 and 3π/16 reminiscent of, but distinctly different from, the Hofstadter-like spectral butterflies observed in the spectra for twisted bilayers and Bloch electrons in magnetic fields 21 . The distinct anomaly in the other figure panels associated with this "butterfly" is due to a region of parameter space associated with very short system period. A careful comparison of the visual features between the eigenvalue density and eigenvector IPR strongly suggests significant correlations between the eigenvalues and eigenvectors.
The eigenvector expansion of χ 1 E in Eq. (6) provides a clear connection between resonant frequencies and large field intensity when Im sðωÞ ( 1. However, our analysis of Fig. 3 also indicates these resonant frequencies correspond to fields that are either extended throughout the medium, as in the leftmost panel of Fig. 3a, or to a mixture of localized and extended states giving rise to more spatially varied field characteristics in both the intensity and localization, as in the leftmost panel of Fig. 2b, with sensitive dependence on frequency. We now make this correspondence more precise in an analysis of the magnitude and phase of ϵ * and the localization of E and D. They are displayed in Fig. 4 for various values of the Moiré twist angle θ, for 0°≤ ϕ ≤ 2°, as a function of Re sðωÞ. The Drude peak at ω = 0 (s(0) = 0) present for all values of ϕ indicates the composite is conducting for ω = 0 39 . For ϕ = 0, at the resonant frequencies both μ and |ϵ * | are sharply peaked and ϵ * diverges as Im s ! 0. These frequencies correspond to the so-called surface plasmon resonance, which characteristically shows up as a strong absorption line in experiments 39 . At these resonant frequencies ϵ * also undergoes a dramatic switch in phase which gives rise to an "optical transition," where the material response changes from inductive (metallic) to capacitive (dielectric)a phenomenon observed in optical cermets 40 . These phase switches also occur at the troughs of |ϵ * |, where |ϵ * | and the mass of μ are small. At these band gap frequencies the material behaves effectively like an electrical insulator. As ϕ increases, the transition frequencies still correspond to the peaks and troughs in |ϵ * |, though the frequency dependence of these features becomes more irregular.
The IPR for |χ 1 E| (normalized to unit length) provides a measurement of localization for the electric field itselfequivalently for the normalized displacement field χ 1 D = ϵ 1 χ 1 E.  ) and (e), slightly saturated at the ends to reveal more detail). Short period systems appear as horizontal streaks; for ρ and μ sharp isolated resonances are identified by localized yellow resonant peaks surrounded by dark blue troughs with values orders of magnitude smaller. The influence of μ on ϵ * /ϵ 2 is clear, with striking similarities. For the IPR in (e) extended and localized vectors are identified by dark blue (with GOE value labeled) and yellow, with mobility edges indicated by sudden changes from one extreme to the other. Some of the θ values associated with these short period systems are identified by black tick marks on the right, labeled by the bound K ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 þ n 2 p on the system period. Figures 3c and 4b show there is a close relationship between the eigenvector IPR, IPR(v), plotted versus λ j and the electric field IPR, IPR(E), plotted versus Re sðωÞ, as anticipated above. Specifically, there are frequency regions where the eigenmodes and the electric field are simultaneously localized or extended. Moreover, for ϕ = 0, 1/8, and 1/2 there are several clear mobility edges in IPR(E), following those in IPR(v), showing high variability in field localization for small changes in s(ω), which also correspond to resonant frequencies and high variability in field intensity.
Physical implications. In Figure 2b the localized (loc) and extended (ext) fields for ϕ = 1/8, 1/2, and 2 were computed for values of Re sðωÞ with optical frequencies ωindicated by red dots in Figure 4. Comparing these two figures for the panels with values ϕ = 1/8 and 1/2 further demonstrates the frequencydependent localization/delocalization transition in the displacement field for the same microstructure. Moreover, the panels for localized (loc) fields in Figure 2b also correspond to resonant peaks in μ in Fig. 3b, which accounts for the high variability in the field intensity in Figure 2b and the amplitude of ϵ * in Fig. 4a. Furthermore, Fig. 4 for ϕ = 1/8 and 1/2, shows that toward the infrared end of the spectrum the displacement field is extended and the response of ϵ * is inductive (metallic), while toward the ultraviolet end of the spectrum the displacement field is more localized and the response of ϵ * is capacitive (dielectric). There are also band gap frequencies in the optical range. As ϕ surpasses 1/8, band gap frequencies are absent. The larger checkerboard scale for |χ 1 E| shown in Fig. 2b decreases in size and all the material characteristics described above begin to qualitatively resemble those of the random percolation model for p = p c as ϕ → 2. The more regularly distributed eigenvector localization gives rise to spatially varied, meandering, tenuously connected field paths as shown in the corresponding panels of Fig. 3a.
These observations indicate a high degree of tunability in the frequency dependence of the phase and magnitude of ϵ * and the localization and intensity of E and D. The resonant and band gap frequencies present for small ϕ are tunable through the microstructure itself via the scale r and Moiré twist angle θ in Eq. (1). We predict that these material characteristics can be reproduced experimentally and tuned by fabrication methods used for etched metallic substrates. (In ref. 46 , a small change in Moiré twist angle for bilayer graphene induces a change in conductivity similar to what we observe here for ϵ * ). Since the transformation in Eq. (1) is deterministic, one can also obtain material characteristics similar to those of random systems in a predictable, reproducible manner. This tunability makes our Moiré-type composite class an ideal test bed for potential engineering applications.
Random matrix theory analysis. Statistical quantities for the eigenvalues λ j of μ provide insights into why the high-density resonances of μ, present for the short period system with ϕ = 0, spread out as ϕ increases and the system becomes quasiperiodic. The nearest neighbor eigenvalue spacing distribution (ESD) P(z) was initially introduced in random matrix theory to describe fluctuations of characteristic quantities for random systems, but has since accurately described quantities for non-random systems with sufficient complexity 47 . The ESD probes short-range correlations of eigenvalues 47 . For highly correlated Wigner-Dyson (WD) spectra exhibited by, for example, the Gaussian orthogonal ensemble (GOE) of real-symmetric random matrices, the ESD is accurately approximated by PðzÞ % ðπz=2Þ expðÀπz 2 =2Þ, Wigner's surmise, which illustrates eigenvalue repulsion, vanishing linearly as spacings z → 0 47,48 . In contrast, the ESD for uncorrelated Poisson spectra, PðzÞ ¼ expðÀzÞ, allows for significant level degeneracy 47 . Figure 6 a displays the ESD for the eigenvalues λ j of G for several values of 0°≤ ϕ ≤ 2°. The blue dash-dot curve is the ESD for Poisson spectra, while the green dashed curve is the ESD for the GOE. For ϕ = 0, 1/64, and 1/32, the sharply peaked resonances in μ with high eigenvalue density give rise to a significant probability of zero spacings, with P(0) ≳ 0.4. However, as ϕ increases and the composite microgeometry becomes quasiperiodic, the behavior of the ESDs starts to be characterized by weakly correlated Poisson-like statistics 48 , also observed for eigenvalues of G for the low volume fraction percolation model 27 . They increase linearly from zero but the initial slope of P(z) is steeper than in the WD case, implying less level repulsion. As ϕ → 2, the slope of P(z) decreases, indicating an increase in level repulsion, causing the eigenvalues of μ to spread out as the ESD transitions toward obeying that of the GOE, characterized by highly correlated eigenvalues with strong level repulsion.
A broader overview of the Moiré parameter space. We conclude this section with a discussion of Fig. 6b, which displays the average eigenvector IPR with yellow hues corresponding to short period systems with highly extended eigenmodeshence displacement fieldsand mobility edges, and dark green to blue hues corresponding to quasiperiodic, random-like systems with more regularly distributed eigenmodes and meandering, tenuously connected field paths. Our results here are only a snapshot, which nevertheless reveals the great diversity of this class of composite materials with myriad microgeometric variations, each with a potentially distinct frequency dependence in both the phase and magnitude of ϵ * and the localization and intensity of E and D. Figure 1 shows that the arrangement of finite period systems is fractal in nature. It is clear from Figs. 1 and 6b that we have merely scratched the surface in describing this fascinating class of composite materials with tuneable capabilities in both frequency and geometry, potentially enabling materials to be fabricated that achieve desired field characteristics and dielectric responses suitable for a broad range of engineering applications.

Conclusion
A class of Moiré-structured 2D composite materials is introduced. Bulk transport is explored using a Stieltjes integral representation for the effective transport coefficients, and the complex permittivity ϵ * in particular. The representation involves a spectral measure μ of a real-symmetric matrix G, and a summation formula for the displacement field D, involving the eigenvalues λ i and eigenvectors v i of G. The localization properties of D and the dielectric profile for ϵ * are analyzed as the Moiré twist angle θ varies 2 degrees. This small change in θ gives rise to a sharp transition in the microgeometry of the composite material, from periodic to quasiperiodic as the period increases ad infinitum. Short period systems are characterized by sharp resonances in μ which give rise to optical frequencies ω where ϵ * is sharply peaked (so-called surface plasmon resonance frequencies) and ϵ * undergoes an "optical transition" from inductive (metallic) to capacitive (dielectric). Band gap optical frequencies are also observed. Moreover, D is highly extended for certain ranges of frequency, separated by small "mobility edge" frequency regions of large localization variability, that follow the resonant peaks of μ with high intensity regions of D. These characteristics make the dielectric profile and field response highly tunable, a desired feature in engineering applications. As the system is tuned to quasiperiodicity, an increase in eigenvalue repulsion, as measured by the eigenvalue spacing distribution (ESD), causes the sharp resonances of μ to spread out, while the localization characteristics of D and the dielectric profile of ϵ * begin to qualitatively resemble those of the percolation model near its transition point. It is suggested that these material characteristics could be reproduced experimentally and tuned by fabrication methods used for etched metallic substrates.

Code availability
Mathematical and numerical methods used to compute the spectral measures and associated spectral statistics displayed in this manuscript are detailed in ref. 26 . Associated code will be made available upon reasonable request.

Data availability
Numerical data used to generate figures in this manuscript will be made available upon reasonable request.