Introduction

In recent years, the engineering of self-similar structures1 in photonics and nano-optics technologies2,3,4,5,6,7 enabled the manipulation of light states beyond periodic8 or disordered systems9,10, adding novel functionalities to complex optical media11,12,13,14 with applications to nano-devices and metamaterials15,16,17. The concept of multifractality (MF), which describes intertwined sets of fractals18,19,20, was first introduced in order to analyze multi-scale energy dissipation in turbulent fluids20 and broadened our understanding of complex structures that appear in various fields of science and engineering21. Specifically, multifractal concepts provided a number of significant insights into signal analysis22, finance23,24,25, network traffic26,27,28, photonics12,29,30,31,32,33,34,35 and critical phenomena36,37,38,39,40,41,42. In fact, critical phenomena in disordered quantum systems have been the subject of intense theoretical and experimental research leading to the discovery of multifractality in electronic wave functions at the metal-insulator Anderson transition for conductors36,40, superconductors38, as well as atomic matter waves43. The multifractality of classical waves has also been observed in the propagation of surface acoustic waves on quasi-periodically corrugated structures35 and in ultrasound waves through random scattering media close to the Anderson localization threshold44.

Considering the fundamental analogy between the behavior of electronic and optical waves10, the question naturally arises on the possibility to experimentally observe and characterize multifractal optical resonances in the visible spectrum using engineered photonic media. Besides the fundamental interest, multifractal optical waves offer a novel mechanism to transport and resonantly localize photons at multiple-length scales over extended surfaces, enhancing light–matter interactions across broad frequency spectra. These are important attributes for the development of more efficient light sources, optical sensors and nonlinear optical components12,45. However, to the best of our knowledge, the direct experimental observation of multifractal optical resonances in engineered scattering media is still missing.

In this paper, we demonstrate and systematically characterize the multifractal behavior of optical resonances in aperiodic arrays of nanoparticles with the distinctive aperiodic order that is intrinsic to fundamental structures of algebraic number theory46,47,48. The paper is organized as follows. First, we analyze the far-field diffraction spectra of these aperiodic photonic systems, which have been conjectured to exhibit singular-continuous characteristics49. Second, we perform dark-field microscopy measurements to directly visualize the scattering resonances of the fabricated structures across the visible spectrum where they exhibit complex spatial distributions with intensity oscillations at multiple-length scales. Third, we perform frequency–frequency correlation analysis that subdivides the measured scattering resonances according to their common structural features, enabling a classification that greatly reduces the dimensionality of the measured datasets. Finally, we demonstrate that the scattering resonances of the proposed photonic arrays exhibit strong multifractal behavior across the entire visible range. Our findings show the coexistence of resonant states with different localization properties and multifractal intensity fluctuations in deterministic two-dimensional structures based on algebraic number theory. We believe that the results presented here provide yet-unexplored possibilities to exploit multifractality as a novel engineering strategy for optical encoding, sensing, lasing and multispectral devices.

Results and discussion

The investigated devices consist of TiO2 nanopillars deposited atop a transparent SiO2 substrate and arranged according to the prime elements of the Eisenstein and Gaussian integers47, as well as two-dimensional cross sections of the irreducible elements of the Hurwitz and Lifschitz quaternions48. These arrays have recently been shown to support spatially complex resonances with critical behavior49 akin to localized modes near the Anderson transition in random systems40,44. Figure 1a outlines the process flow utilized for the fabrication of the arrays (details can be found in the Methods section). Scanning electron microscope (SEM) images of the fabricated devices are reported in Figs. 1b, e, whereas their geometrical properties, illustrated in Supplementary Figs. 1 and 2, are discussed in the Supplementary Note 1. The multifractality in the optical response of these novel photonic structures is demonstrated experimentally by analyzing the intensity fluctuations of their scattering resonances that are measured using dark-field scattering microscopy across the visible spectral range. Before addressing the multifractal nature of the measured scattering resonances, in the next section we discuss the far-field diffraction properties of the fabricated arrays because they provide insights on the fundamental nature of these novel aperiodic systems.

Fig. 1: Experimental realization of sub-wavelength prime arrays and their diffractive properties.
figure 1

a Process flow used to fabricate the TiO2 nanocylinders on a SiO2 substrate. b–e Scanning electron microscope images of fabricated samples: b–e Eisenstein, Gaussian, Hurwitz and Lifschitz prime configurations, respectively. Insets: enlarged view of the central features for each pattern type. Nanocylinders have a 210 nm mean diameter, 250 nm height and average inter-particle separation of 450 nm. Calculated (f–i) and measured (j–m) k-space intensity profiles corresponding to the prime arrays reported in b–e, respectively. n Optical setup for measuring far-field diffraction patterns of laser illuminated sample. IP image plane, AS Aperture stop, FP Fourier plane.

Full size image

Spectral characterization of fabricated arrays

Aperiodic systems can be rigorously classified according to the nature of their Fourier spectral properties. In fact, according to the Lebesque decomposition theorem50, any positive spectral measure (here identified with the far-field diffraction intensity) can be uniquely expressed in terms of three primary spectral components: a pure-point component that exhibits discrete and sharp (Bragg) peaks51, an absolutely-continuous one characterized by a continuous and differentiable function52 and a more complex singular-continuous component where scattering peaks appear to cluster into a hierarchy of self-similar structures forming a highly-structured spectral background53,54. Systems with singular-continuous spectra do not occur in nature. Moreover, their spectral properties, which are often associated to the presence of multifractal phenomena, are difficult to characterize mathematically51. In this section, we introduce and characterize the spectral properties of our scattering systems by investigating their measured far-field diffraction spectra.

When laser light uniformly illuminates the arrays at normal incidence, sub-wavelength pillars behave as dipolar scattering elements and the resulting far-field diffraction pattern is proportional to the structure factor defined by:

$$S({\boldsymbol{k}}) = \frac{1}{N}\left| {\mathop {\sum}\limits_{j = 1}^N {e^{ - i{\boldsymbol{k}} \cdot {\boldsymbol{r}}_j}} } \right|^2,$$
(1)

where k is the in-plane component of the wavevector and rj are the vector positions of the N nanoparticles in the array. This is demonstrated in Fig. 1, where we compare the calculated structure factors, shown in Figs. 1f–i, with the measured diffraction patterns, shown in Figs. 1j–m. The optical setup used to measure the diffraction is schematically illustrated in Fig. 1n and further discussed in the Methods section. The experimental diffraction patterns match very well with the calculated ones and display sharp diffraction peaks embedded in a weaker diffuse background, particularly noticeable in the case of Gaussian and Eisenstein primes. The coexistence of sharp diffraction peaks with a structured (self-similar) background is characteristic of singular-continuous spectra described by singular functions that oscillate at every length scale49,51,53. The strength of the continuous spectral component weakens progressively from Eisenstein and Gaussian primes to Hurwitz and Lifschitz structures, which display a more regular geometry. See also Supplementary Fig. 3 and the related discussion in the Supplementary Note 1.

A remarkable property of structures that support singular-continuous components is that their resonant modes are critical55,56,57,58. Critical modes exhibit local fluctuations and spatial oscillations over multiple-length scales that are quantitatively described by multifractal analysis45,59. These peculiar scattering resonances have been demonstrated in one-dimensional aperiodic systems55,56,57,58 and in the propagation of surface acoustic waves on quasi-periodically corrugated structures35. In the following, we experimentally demonstrate, using dark-field imaging technique, the multifractal nature of the scattering resonances of aperiodic arrays based on generalized prime numbers in the visible range.

Dark-field scattering measurements

Dark-field (DF) microscopy is a well-established technique to image the intensity distributions of the scattering resonances supported by dielectric and metallic nanoparticles arrays60,61, to visualize the dynamic of microtubules-associated protein62 and to analyze the formation of colorimetric fingerprints of different aperiodic nano-patterned surfaces63,64. In this section, we use DF microscopy to directly visualize the distinctive scattering resonances supported by the fabricated arrays of nanoparticles at multiple wavelengths across the visible spectrum. In particular, we use a confocal DF microscopy setup in combination with a line-scan- hyperspectral imaging system in illumination-detection configuration as explained in the Methods section (see also Supplementary Note 2). A large dataset of 355 different dark-field scattering images, for each investigated structure, was collected corresponding to the excitation of the scattering resonances at different frequencies ranging from 333 THz up to 665 THz with a resolution of ~0.9 THz.

Figure 2 shows representative dark-field images of the scattering resonances of Eisenstein  2a–d, Gaussian  2e–h, Hurwitz  2i–l and Lifschitz  2m–p arrays collected at different frequencies across the visible range, as specified by the markers in Fig. 3. The scattering resonances of the Eisenstein and Gaussian prime structures display a clear transition from localization in the center of the arrays to a more extended nature in the plane of the arrays. On the other hand, this scenario is far richer for the Hurwitz and Lifschitz configurations. In particular, we found that the spatial distributions of their scattering resonances exhibit the following characteristics: (i) weakly localized around the central region of the arrays, (ii) localized at the edges of the arrays and (iii) spatially extended over the whole arrays (more details are provided in Supplementary Note 3, Supplementary Fig. 4 and in the next section).

Fig. 2: Direct observation of the scattering resonances of prime arrays.
figure 2

Representative multispectral dark-field images of the scattering resonances of Eisenstein (a–d), Gaussian (e–h), Hurwitz (i–l) and Lifschitz (m–p) prime arrays. These quasi-modes are the scattering resonances that most interact with their corresponding geometrical supports at a fixed frequency. The scale bars indicate the scattering resonance intensities normalized with respect to the transmission signal of a reference Ag-mirror and reported in arbitrary units. The x and y axis express the length and the width of the investigated devices expressed in micron scale.

Full size image
Fig. 3: Frequency–frequency correlation analysis of scattering resonances.
figure 3

a–d Frequency–frequency correlation matrix C(ν, ν′) of Eisenstein, Gaussian, Hurwitz and Lifschitz arrays, respectively. The white line indicates the diagonal elements, i.e., the normalized variance C(ν, ν) of the scattered intensity. The different markers identify the frequency location of representative scattering resonances that capture the main features of the correlation analysis. e–h Normalized variance as compared to the normalized density of states (DOS) as a function of the frequency \(\nu\).

Full size image

All the measured scattering resonances exhibit complex spatial distributions characterized by intensity fluctuations at multiple-length scales. However, by inspecting Fig. 2, we can also recognize structural features in the intensity distributions that are shared by modes of different frequencies (compare for instance Fig. 2j with 2l). In order to quantitatively classify the measured scattering resonances according to their common structural features, we employ in the next section frequency–frequency correlation analysis.

Correlation analysis of dark-field scattering resonances

Image cross-correlation techniques are widespread in signal processing where they are used to compare different data and identify structural similarities. In this section, we apply the frequency–frequency correlation analysis introduced by Riboli et al.65 to a dataset of 1420 measured scattering resonances. This approach enables to reduce the enormous complexity of our total dataset into a few constituent groups of resonances that are identified by the frequency–frequency correlation matrix. The off-diagonal elements of this matrix provide the degree of spatial similarity between two scattering resonances with frequencies \(\nu\) and \(\nu\)′, as defined below65:

$$C(\nu ,\nu ^\prime ) = \frac{{\frac{1}{n}\mathop {\sum}\limits_{s = 1}^n {\delta I} ({\boldsymbol{r}}_s,\nu )\delta I({\boldsymbol{r}}_s,\nu ^\prime )}}{{\left[ {\frac{1}{n}\mathop {\sum}\limits_{s = 1}^n I ({\boldsymbol{r}}_s,\nu )} \right]\left[ {\frac{1}{n}\mathop {\sum}\limits_{s = 1}^n I ({\boldsymbol{r}}_s,\nu ^\prime )} \right]}},$$
(2)

where the counting variable s identifies the pixel of the intensity I of the considered dark-field image, n is the total number of pixels, rs = (xs, ys) denotes the in-plane coordinates, and \(\delta I({\boldsymbol{r}}_s,\nu ) = I({\boldsymbol{r}}_s,\nu ) - [\mathop {\sum}\nolimits_{s = 1}^n I ({\boldsymbol{r}}_s,\nu )]{\mathrm{/}}n\) describes the fluctuations of the intensity I (see also the Methods section and Supplementary Note 3).

The results of the correlation analysis separate the set of collected scattering resonances into two spectral classes. The first class includes the Eisenstein and Gaussian arrays and is characterized by large fluctuations in the correlation matrix concentrated around two distinct spectral regions, as shown in Figs. 3a, b. Specifically, these two well-defined spectral regions are located at \(\nu\) ≈ 460 THz and \(\nu\) ≈ 600 THz, respectively. The off-diagonal peaks of the correlation matrices reported in Figs. 3a, b show that these resonances are correlated, indicating that their spatial distributions share common structural features. Moreover, a reduction of the intensity fluctuations is observed around \(\nu\) ≈ 550 THz in both Eisenstein and Gaussian arrays. All these features can also be visually recognized in the spatial distributions of the selected scattering resonances reported in Figs. 2a–h. In fact, the representative resonances displayed in Fig. 2, which correspond to the different markers in Figs. 3a–d, have been identified based on the results of the correlation analysis. Specifically, Figs. 2a, c do not share any common features: the resonance in Fig. 2a shows the maximum of its intensity at the center of the array, whereas the resonance shown in Fig. 2c displays an intensity minimum in the same region. On the other hand, the spatial distributions displayed in Figs. 2b, d are highly correlated. Similar results are obtained for the Gaussian arrays: while the resonances reported in Fig. 2e and in Fig. 2g are not correlated, the spatial distributions in Figs. 2f, h exhibit a similar structure.

In the second class, which includes the Hurwitz and Liftschitz arrays, the C(\(\nu\), \(\nu\)′) matrices are more structured and exhibit multi-scale fluctuations spreading over the entire frequency spectrum, as shown in Figs. 3c, d. Specifically, their normalized variances are characterized by peaks and dips that spread over the entire measured spectrum. This characteristic behavior indicates that the spatial intensity distributions of the measured scattered radiation rapidly fluctuate with frequency due to the excitation of scattering resonances in the systems. In the Hurwitz configuration, for example, some of these fast fluctuating resonances are correlated (or anti-correlated), as highlighted by the off-diagonal elements of the correlation matrix of Fig. 3c. In fact, the scattering resonances of the Hurwitz array can be classified in four different spectral regions: (i) [350, 420] THZ, (ii) [430, 480] THz, (iii) [480, 600] THz, (iv) larger than 600 THz (see also Supplementary Fig. 5). Specifically, region (i) is characterized by the same scattering resonances of Fig. 2i that are mostly localized at the center of the array and that are anti-correlated with all the others. On the other hand, modes in region (ii) are correlated with the ones in region (iv). This correlation can also be recognized by visually comparing the scattering resonance of Fig. 2j with the one reported in Fig. 2l. Finally, region (iii) is populated by scattering resonances spatially localized at the edge of the arrays, as shown in Fig. 2k. The same considerations apply to the Lifschitz prime array, as shown in Figs. 2m–p and 3d. This complex scenario is well-reproduced also in our numerical simulations, shown in Supplementary Figs. 6 and 7. The scattering resonances in the second class display a more complex spatial structure as compared to the ones of the first class. This is confirmed by comparing the normalized variance of the frequency–frequency correlation matrices with the spectral behavior of the computed optical density of states (DOS) of the systems, which we obtained using the Green’s matrix spectral method (see the detailed derivation in the Supplementary Note 4).

Figures 3e–h present the comparison between the DOS and the normalized variance C(\(\nu\), \(\nu\)) (see Methods for more details). The Eisenstein and Gaussian prime arrays feature a DOS with two main peaks demonstrating that their scattering resonances are predominantly located within the two spectral regions identified using the correlation analysis. On the other hand, the behavior of C(\(\nu\), \(\nu\)) and of the DOS for the Hurwitz and Lifschitz configurations is far richer and features multiple spectral regions of strong intensity fluctuations. Although this comparison remains qualitative in nature, it is consistent with the presence of spectral sub-structures that appear on a finer scales in these arrays, which is typical of aperiodic systems with singular-continuous spectra51,53,59. Interestingly, the classification of our structures into two broad spectral classes may reflect a fundamental number-theoretic difference in the nature of the corresponding rings. In fact, although the Eisenstein and Gaussian structures are constructed based on the prime elements of the rings associated to the imaginary quadratic fields \({\Bbb Q}(\sqrt { - 3} )\) and \({\Bbb Q}(\sqrt { - 1} )\), which are the commutative rings \({\Bbb Z}[( - 1 + i\sqrt 3 )/2]\) and \({\Bbb Z}[i]\), the Hurwitz and Lifschitz structures are constructed based on two-dimensional cross sections of the irreducible elements of quaternions, which form a non-commutative ring66,67,68.

Multifractal analysis of dark-field scattering resonances

In this section, we directly demonstrate that the measured dark-field scattering resonances of the investigated arrays exhibit strong multifractal behavior. We limit our analysis to the representative scattering resonances shown in Fig. 2 that, as previously discussed based on the frequency–frequency correlation matrix, are representative of the variety of structural features observed across the visible spectrum. In order to quantitatively describe intensity oscillations that occur at multiple scales, we apply the multifractal scaling method developed by Chhabra et al.19. In particular, we employ the box-counting method to characterize the size-scaling of the moments of the light intensity distribution. This is achieved by dividing the system into small boxes of varying size l. We then determine the minimum number of boxes N(l) needed to cover the system for each size l and evaluate the fractal dimension Df using the power-law scaling \(N(l) \sim l^{ - D_f}\).

Traditional (homogeneous) fractal structures are characterized by a global scale-invariance symmetry described by a single fractal dimension Df. On the other hand, heterogeneous fractals or multifractals are characterized by a continuous distribution \(f(\alpha)\) of local scaling exponents α such that \(N(l) \sim l^{ - f(\alpha )}\)19. The so-called singularity spectrum f(α) generalizes the fractal description of complex systems in terms of intertwined sets of traditional fractal objects. Different approaches are used to extract f(α) from the local scaling analysis.

Here, we have obtained the multifractal spectra directly from the dark-field measurements by employing the method introduced in ref. 19. Details on the implementation are discussed in the Supplementary Note 5 and in Supplementary Figs. 8–11. Figures 4a–d demonstrate the multifractal nature of the measured scattering resonances of the prime arrays. All the singularity spectra f(α) exhibit a downward concavity with a large width Δα, which is the hallmark of multifractality19. The behavior of the width Δα as a function of frequency reflects the previously identified classification in terms of the correlation matrix (Supplementary Fig. 10). In particular, the arrays with more significant intensity fluctuations, i.e., Eisenstein and Gaussian arrays, also display the broadest singularity spectra. Therefore, our findings establish a direct connection between the multifractal properties of the geometrical supports of the arrays (discussed in the Supplementary Note 1), the MF of the corresponding scattering resonances and the singular-continuous nature of their diffraction patterns. To further characterize the MF of the scattering resonances, we also compute, based on the experimental data, the mass exponent τ(q) and the generalized dimension D(q), which provide alternative descriptions of multifractals. Multifractals are unambiguously characterized by a nonlinear τ(q) function and a smooth D(q) function18. This is reported in Supplementary Fig. 9 and discussed in the Supplementary Note 6.

Fig. 4: Multifractal nature of scattering resonances of prime arrays.
figure 4

Multifractal singularity spectra f(α) of the local scaling exponent α (a–d), and probability density functions \(P({\mathrm{ln}}\,\bar I_{\mathrm{l}})\) of the logarithm of the box-integrated intensity \(\bar I_{\mathrm{l}}\) (e–h) of representative scattering resonances of Eisenstein (a, e), Gaussian (b, f), Hurwitz (c, g) and Lifschitz (d, h) prime arrays, respectively. Markers represent the data and the continuous lines in e–h show the best fits obtained using log-normal function model. Different colors refer to different frequencies corresponding to the main features of their correlation matrices. Error bars in a–d represent the standard deviations.

Full size image

A direct consequence of multifractality is the non-Gaussian nature of the probability density function (PDF) of the light intensity distribution near, or at, criticality36,44. We have evaluated the PDF of the scattered radiation from the histogram of the logarithm of the box-integrated intensities, i.e., \(P({\mathrm{ln}}\,\bar I_{\mathrm{l}})\). We show in Figs. 4e–h the histograms produced by a box size of ~0.8 × c/\(\nu\), where c is the speed of light and \(\nu\) is the frequency of the considered scattering resonance. We note that, although the intensity distributions of the scattering resonances in Fig. 2 show different degrees of spatial variations, all the PDFs are very-well reproduced by a log-normal function model (continuous lines in Figs. 4e–h). Furthermore, we have verified that these findings do not depend on the box-counting size l, as shown in Supplementary Fig. 11.

Finally, the multifractal spectrum f(α) can be rigorously obtained from the PDF of the intensity within the parabolic approximation21,44:

$$f(\alpha ) = D_f - (\alpha - \alpha _0)^2{\mathrm{/}}4(\alpha _0 - d),$$
(3)

where Df = f(α0) is the fractal dimension of the geometrical support and d is the dimensionality of the investigated problem (d is equal to 2 in our analysis). Equation (3) uniquely associates f(α) to the PDF of the logarithm of the box-integrated intensities through the parameter α0 via the formula \(\alpha _0 = 2\mu d{\mathrm{/}}(2\mu + \sigma ^2)\), where μ and σ are the mean and the standard deviation of the PDF, respectively. In our data, we find deviations from the parabolic spectrum, which are associated to the strong MF of the scattering resonances of prime arrays across a broad range of frequencies44 (see Supplementary Fig. 12 and the related discussion in Supplementary Note 7). The observed strong MF in the scattering resonances is thus uniquely associated to the multi-scale geometrical structure of prime arrays. In contrast, only weak MF was previously reported in uniform random scattering media at their metal-insulator Anderson transitions36,44.

In conclusion, we have provided an experimental observation of multifractality of classical waves in the optical regime by studying the scattering resonances of deterministic photonic structures designed to reflect the intrinsic aperiodic order of complex primes and irreducible elements of quaternion rings. Beyond its fundamental interest with respect to the discovery and characterization of multifractality in certain fundamental structures of number theory69,70, our findings establish a novel mechanism to transport and resonantly localize photons at multiple-length scales and to enhance light–matter interactions, with applications to active photonic devices and novel broadband nonlinear components. Moreover, the concepts developed in this work can be naturally transferred to quantum waves71,72 and stimulate the exploration of novel quantum phases of matter and many-body phenomena that emerge from fundamental structures of algebraic number theory.

Methods

Sample fabrication

TiO2 thin films were grown by reactive direct current (DC) magnetron sputtering (MSP) on quartz (SiO2) substrates with a Denton Discovery D8 Sputtering System. A 3 inch diameter Ti target (Kurt J. Lesker, 99.998%) was used. The deposition was performed under a 200 Watt power at a 3.0 mTorr deposition pressure with a 1:3 Sccm flow rate ratio of O2 to Ar gases. These growth conditions resulted in a deposition rate of about 1 nm/min. Chamber base pressure was kept below 5 × 10−7 Torr, the target–substrate distance was fixed at 10 cm, and substrates were rotated at a speed of 0.1 Hz. Substrates were solvent washed and plasma cleaned in O2 prior to deposition.

The optical constants of the films were determined using a variable angle spectroscopic ellipsometer (V-Vase, J.A. Woollam) in the wavelength range of 300–2000 nm. The measured data were fitted in good agreement with the Cauchy model, which provides an empirical relationship between refractive index and wavelength. The film thickness, verified with ellipsometry and SEM, was targeted at 250 nm.

Nanoparticle fabrication was performed with electron beam lithography. A polymethyl methacrylate (PMMA) resist layer of about 100 nm thickness was spun and baked before sputtering a thin conducting layer (~6 nm) of Au to compensate for the electrically insulating SiO2 substrate. The resist was exposed at 30 keV by an SEM (Zeiss Supra 40) integrated with Nanometer Pattern Generation System direct write software. The sample was then developed in isopropyl alcohol to methyl isobutyl ketone (IPA:MIBK) (3:1) solution for 70 s followed by a rinsing with IPA for 20 s. Next, a 20-nm-thick layer of Cr was deposited on top of the developed resist with electron beam evaporation (CHA Industries Solution System). After this deposition, unwanted Cr was removed with a three minutes lift-off in acetone and the nanoparticle patterns were transferred from the Cr mask to the TiO2 thin film by reactive ion etching (RIE, Plasma-Therm, model 790) using Ar and SF6 gases. Finally, the Cr mask was removed by wet etching in Transene 102073.

Far-field diffraction setup

We measure the far-field diffraction pattern of laser light passing through the sample using the setup depicted in Fig. 1n. A 405 nm laser is focused onto the patterned area to overfill the pattern and to ensure uniformity in intensity across it. The forward scattered light is collected by a high numerical aperture objective (NA = 0.9 Olympus MPlanFL N), which collects light scattered up to 64° from the normal direction. A 4-F optical system, immediately behind the objective, creates an intermediate image plane and intermediate Fourier plane. An iris located at the intermediate image plane was used to restrict the light collection area only to the patterned region. A second 4-F optical system then re-images the intermediate Fourier plane onto the charged coupled device (CCD) with the appropriate magnification. Finally, digital filtering was employed to remove the strong direct component of the diffraction spectra to produce clear images.

Dark-field scattering setup and image acquisition

In order to measure the spatial distribution of the scattering resonances of prime arrays, we used a line-scan-hyperspectral imaging system (XploRA Plus by Horiba Scientific). This integrated microscope platform provides both widefield and hyperspectral imaging. Its enhanced dark-field system produces much higher signal-to-noise than standard dark-field imaging systems, such the one described in Supplementary Note 2, which results in much cleaner and accurate images, optimizing dark-field detection capability for non-fluorescing nanoscale sample. Specifically, the field of view is illuminated by light from a lamp through the microscope objective. A line-like region of the field of view is projected into the spectrograph entrance slit. The slit image is then spectrally dispersed by means of a grating onto the camera image sensor. This allows the collection of all the spectra corresponding to a single line in the field of view. Scanning the sample with a stage then produces the hyperspectral image of the entire field of view. In our system, scanning 696 lines results in square-size hyperspectral image. The resolution of the CCD (Cooke/PCO Pixelfly PCI Camera) was set to be 430 × 470, where 430 and 470 corresponds to a spatial and spectral resolution of 30 nm and 0.25 THz, respectively. The exposure time was set to be 1000 ms with a time interval of 1 min between different scanning lines. Finally, the dark-field data were normalized with respect to a reference signal of an Ag-mirror.

Correlation analysis

Each element C(\(\nu\)i, \(\nu\)j) of the frequency–frequency correlation matrices of Figs. 4a–d is the result of an average over 9 × 104 correlated values. Equation (2) can be written in the compact form65:

$$C(\nu ,\nu ^\prime ) = \frac{{\langle \delta I({\boldsymbol{r}},\nu )\delta I({\boldsymbol{r}},\nu ^\prime )\rangle }}{{\langle I({\boldsymbol{r}},\nu )\rangle \langle I({\boldsymbol{r}},\nu ^\prime )\rangle }}.$$
(4)

The normalization of the covariance \(\langle \delta I({\boldsymbol{r}},\nu )\delta I({\boldsymbol{r}},\nu ^\prime )\rangle\) with respect to the product of the average values \(\langle I({\boldsymbol{r}},\nu )\rangle \langle I({\boldsymbol{r}},\nu ^\prime )\rangle\) minimizes the intrinsic spectral effects related to the illumination source of the experimental apparatus. Each element of \(C(\nu, {\nu}^{\prime})\) is a combination of spatial intrinsic fluctuations of the system’s parameters, extrinsic effects related to the illumination-collection efficiency, and to the point spread function of the experimental apparatus. The spectral dependence of the extrinsic effects can be mitigated by the proper normalization of the correlations matrix65. The white diagonal elements of Figs. 3a–d as well as the continuous lines in Figs. 3e–h are the normalized variance of the scattering resonance intensities of prime arrays and are calculated by the equation:

$$C(\nu ,\nu ) = \frac{{\langle \delta I({\boldsymbol{r}},\nu )^2\rangle }}{{\langle I({\boldsymbol{r}},\nu )\rangle ^2}} = \frac{{\sigma ^2(\nu )}}{{\mu ^2(\nu )}},$$
(5)

where μ(\(\nu\)) is the average value and σ(\(\nu\)) is the standard deviation.

Multifractal analysis

The multifractal analysis of both structural and dynamical properties was performed from the corresponding 600 dpi bitmap image using the direct Chhabra–Jensen algorithm19 implemented in the routine FracLac (ver. September 2015) developed for the National Institutes of Health (NIH) distributed Image-J software package74. This method is a useful tool to determine the scaling properties of a certain image, but has to be handled with care. First of all, the size of the boxes must satisfy some requirements. In the FracLac routine, the largest box should be larger than 50% but not exceed the entire image, whereas the smallest box is chosen to be the point at which the slope starts to deviate from the linear regime in the log(N) versus log(1/r) plot74. Furthermore, the scattering resonance maps are not binary. Therefore, it is necessary to specify the threshold value above which the pixels are part of the object under analysis. Different calculations were performed for several threshold percentages (between 55% and 75%) of the maximum intensity of each scattering resonance of Fig. 2. Another source of error could be the scaling method used during the analysis. For each dark-field image, we employed a linear, a rational and a power scaled series of box sizes74. All these aspects do not affect the main results of our paper, as shown by the error bars of Figs. 4a–d.