## Order to disorder in quasiperiodic composites

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 *R**L**C* 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=\sqrt{10}/3\) and \(\theta =\arctan (1/3)+\phi\) 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 microgeometry — from 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\le {{{{{{{\rm{Re}}}}}}}}\,s(\omega )\le 1\) with \({{{{{{{\rm{Im}}}}}}}}\,s(\omega )\ll 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 < {{{{{{{\rm{Re}}}}}}}}\,s\le 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 *R**L**C* 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 \({\omega }_{p}^{2}=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 *R**L**C* 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 *R**L* − *C* network^{40,41}. The dependence of *s*(*ω*) on frequency *ω* is model specific. For the *R* − *C* and *R**L* − *C* models, \(0\le {{{{{{{\rm{Re}}}}}}}}\,s \, < \,1\) for 0 ≤ *ω* < *∞* and \(\delta \le {{{{{{{\rm{Im}}}}}}}}\,s \, < \,0\), where ∣*δ*∣ can be chosen as small as desired, with \({{{{{{{\rm{Re}}}}}}}}\,s\to 1\) and \({{{{{{{\rm{Im}}}}}}}}\,s\to 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\le {{{{{{{\rm{Re}}}}}}}}\,s\le 1\) and \({{{{{{{\rm{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 \({{{{{{{\rm{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 \({{{{{{{\rm{Im}}}}}}}}\,s(\omega )\ll 1\) in mind, we call resonant frequencies the values of *ω* where \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\approx {\lambda }_{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.

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 *R**L* − *C* model^{40,41}. The dielectric resonances observed for the *R**L* − *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 \({{{{{{{\rm{IPR}}}}}}}}({{{{{{{\bf{u}}}}}}}})={\sum }_{i}{u}_{i}^{4}\), where *u*_{i} is the *i*th 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 \, < \,{{{{{{{\rm{Re}}}}}}}}\,s(\omega )\le 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 applications — given 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=\sqrt{5}/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 density — also 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 \({{{{{{{\rm{Re}}}}}}}}\,s\approx 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 \({{{{{{{\rm{Im}}}}}}}}\,s(\omega )\ll 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 \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\). 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 \({{{{{{{\rm{Im}}}}}}}}\,s\to 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 itself – equivalently for the normalized displacement field *χ*_{1}**D** = *ϵ*_{1}*χ*_{1}**E**. 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 \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\), 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 \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\) 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 frequency-dependent 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)\approx (\pi z/2)\exp (-\pi {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 eigenmodes — hence displacement fields — and 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.