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 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., ∣χ1E∣, calculated via Eq. (6). Since χ1D = ϵ1χ1E 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.
Moiré interference patterns generated by the transformation Eq. (1) give rise to a large class of composite materials with periodic and quasiperiodic microgeometry. (a) Various example subsections arising from different Moiré parameter pairs (r, θ). Microgeometries in (a) are shown with square system sizes 53 for the far left and 73 for all others, small enough to resolve the small-scale detail yet illustrate the large geometric variety. Geometry is illustrated by rendering χ1 = 1 bonds and omitting χ1 = 0 bonds. Cool and warm colors correspond to near-zero and large values of electric field ∣χ1E∣ or displacement ∣χ1D∣, respectively, with the color bar at the top showing the saturated linear scale, normalized to the unit interval. (b) Anderson localization of fields in quasiperiodic media. Departure from short period system geometry is parameterized by ϕ. Composite microgeometry parameterized by \(r=\sqrt{10}/3\) and \(\theta =\arctan (1/3)+\phi\) 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 ∣χ1E∣ 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 = pc27). The inverse participation ratio (IPR) providess quantitative description of this localization phenomenon.
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 = pc.
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 RLC impedance networks40. 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 other39, 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 RLC impedance networks really do resemble the continuum composites they are intended to model39.
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 network40,41. The dependence of s(ω) on frequency ω is model specific. For the R − C and RL − 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 ω increases40,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 interval45.
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 mj 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 mj.
Composite microgeometry and fields, spectral measure μ, and eigenvector inverse participation ration (IPR) plotted for various values of the Moiré twist angle θ, for 0° ≤ ϕ ≤ 2°, scale parameter \(r=\sqrt{10}/3\), and square system of size 73 in (a) to show detail and of size 199 in (b) and (c). The color bars in the upper left of the panels are for reference and indicate \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\) corresponding to impedances throughout the optical frequency range for a Drude model of a gold/vacuum composite. (a) Composite microgeometry and field intensity. Cool and warm colors correspond to near-zero and large values of electric filed ∣χ1E∣ or displacement ∣χ1D∣, respectively, with color bar at the top showing the saturated linear scale. (b) Masses mj of measure μ plotted versus eigenvalues 0≤λj≤1 of the matrix G. Red dots indicate the largest masses, used as indicators in (c) below. (c) IPR values for the eigenvectors vj of G, IPR(v), plotted versus λj. Green and black horizontal lines indicate IPR values for Gaussian orthogonal ensemble (GOE) vectors and completely extended vectors, 3/N1 and 1/N1, respectively. These quantities for the random percolation model at the percolation transition, p = pc = 1/2, are shown in the rightmost panels for comparison.
Relative effective complex permittivity ϵ*/ϵ2 and inverse participation ratio (IPR) of the modulus of the electric field, ∣χ1E∣, normalized to unity, IPR(E), plotted versus \(0\le {{{{{{{\rm{Re}}}}}}}}\,s\le 1\) for \({{{{{{{\rm{Im}}}}}}}}\,s=0.001\) and various values of the Moiré twist angle θ and \(r=\sqrt{10}/3\), for change in twist angle 0° ≤ ϕ ≤ 2°. The color bars in the upper left of the panels are for reference and indicate \({{{{{{{\rm{Re}}}}}}}}\,s(\omega )\) for the optical frequency range for impedances corresponding to the Drude model for a gold/vacuum composite. (a) Amplitude and phase of ϵ*/ϵ2. (b) IPR(E) or equivalently IPR(D). These quantities for the random percolation model at the percolation transition, p = pc = 1/2, are shown in the rightmost panels for comparison. The red dots in (a) and (b) identify values of dependant variable \({{{{{{{\rm{Re}}}}}}}}\,s\) used in Fig. 2b: for ϕ = 1/8, \({{{{{{{\rm{Re}}}}}}}}\,s=0.063,\,0.115\), for ϕ = 1/2, \({{{{{{{\rm{Re}}}}}}}}\,s=0.055,\,0.111\), and for ϕ = 2, \({{{{{{{\rm{Re}}}}}}}}\,s=0.111\).
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 = pc , 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 RL − C model40,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 films41. 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 N1-dimensional unit vector u it is given by \({{{{{{{\rm{IPR}}}}}}}}({{{{{{{\bf{u}}}}}}}})={\sum }_{i}{u}_{i}^{4}\), where ui is the ith component of the vector u, i = 1, …, N1, and satisfies IPR(u) = 1 for a completely localized vector with only one non-zero component and IPR(u) = 1/N1 for a completely extended vector with all components equal in value27. For matrices in the Gaussian orthogonal ensemble (GOE), the eigenvectors are quite extended with a mean asymptotic IPR value of IPRGOE = 3/N127.
Figure 3c displays IPR(vj) for the eigenvectors vj, j = 1, …, N1, 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(vj) have large variability for small changes in λj. As ϕ increases and the microgeometry becomes quasiperiodic, the mobility edges diminish as the values IPR(vj) 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(vj) resemble those of the random percolation model at its threshold p = pc = 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 mj ≳ 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 model26,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 fields21. 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.
(a) Eigenvalue density ρ(λ, θ) (a histogram representation of the density of states ∑jδ(λ − λj)/N1), (b) the spectral function μ(λ, θ) (a kernel estimate representation of the spectral measure), (c) magnitude and (d) phase of relative effective complex permittivity ϵ*/ϵ2, and (e) a histogram-like representation of the inverse participation ratio (IPR) (median IPR of eigenstates associated with each bin – to distinguish mobility edges), all plotted versus the Moiré twist angle θ for scale parameter \(r=\sqrt{5}/2\). We plot these quantities for one full period 0 ≤ θ ≤ π/4. (a), (b), and (e) are plotted versus eigenvalue 0 ≤ λ ≤ 1, while (c) and (d) are plotted versus 0 ≤ Re s≤1 for \({{{{{{{\rm{Im}}}}}}}}\,s=0.001\). Low and high density for ρ and μ and are indicated by dark blue and yellow, respectively, as shown by the color bars (with linear scale in (a), (c), and (e) and log10 scale in (b) 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=\sqrt{{m}^{2}+{n}^{2}}\) on the system period.
The eigenvector expansion of χ1E 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 ω = 039. 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 experiments39. 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 cermets40. 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 ∣χ1E∣ (normalized to unit length) provides a measurement of localization for the electric field itself – equivalently for the normalized displacement field χ1D = ϵ1χ1E. 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 ∣χ1E∣ 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 = pc 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 complexity47. The ESD probes short-range correlations of eigenvalues47. 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 → 047,48. In contrast, the ESD for uncorrelated Poisson spectra, \(P(z)=\exp (-z)\), allows for significant level degeneracy47.
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 statistics48, also observed for eigenvalues of G for the low volume fraction percolation model27. 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) The eigenvalue spacing distribution (ESD) P(z) for various values of (Moiré) twist angle θ, for 0° ≤ ϕ ≤ 2°. The short period system for ϕ = 0 and those with small twist angles 0 ≤ ϕ ≤ 1/32 are characterized by spectral measures μ with very sharp resonances leading to P(0) ≳ 0.4. However, for ϕ ≥ 1/16 the system begins to transition towards obeying Wigner-Dyson statistics with level repulsion, so that P(0) = 0. Level repulsion increases with increasing ϕ as the ESD approaches the Wigner-Dyson ESD, characterized by strong correlations and strong eigenvalue repulsion. (b) The ratio of average eigenvector \(\overline{IPR}\) with IPR GOE = 3/N1 is plotted versus \((r,\tan \theta )\). Yellow hues correspond to short period systems similar to the leftmost panel in Figure 2, characterized by highly extended eigenmodes (hence extended electric and displacement fields) and “mobility edges” with large localization variability. Dark green to blue hues correspond to quasiperiodic systems similar to the one shown in Fig. 2 for ϕ = 2 with material properties that resemble that of random systems with regularly distributed IPR values and tenuously connected electric and displacement field paths. This panel indicates periodic systems have a repeating pattern that turns out to be fractal in nature, as indicated in Fig. 1. Moreover, quite small changes in the Moiré parameters \((r,\tan \theta )\) result in transitions from ordered periodic systems to disordered quasiperiodic, random-like systems.
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.