Multifractality of light in photonic arrays based on algebraic number theory

Many natural patterns and shapes, such as meandering coastlines, clouds, or turbulent flows, exhibit a characteristic complexity that is mathematically described by fractal geometry. Here, we extend the reach of fractal concepts in photonics by experimentally demonstrating multifractality of light in arrays of dielectric nanoparticles that are based on fundamental structures of algebraic number theory. Specifically, we engineered novel deterministic photonic platforms based on the aperiodic distributions of primes and irreducible elements in complex quadratic and quaternions rings. Our findings stimulate fundamental questions on the nature of transport and localization of wave excitations in deterministic media with multi-scale fluctuations beyond what is possible in traditional fractal systems. Moreover, our approach establishes structure–property relationships that can readily be transferred to planar semiconductor electronics and to artificial atomic lattices, enabling the exploration of novel quantum phases and many-body effects. Emergent multifractality is the object of both fundamental and technology-oriented research. Here, the authors demonstrate and characterize multifractality in the optical resonances of aperiodic arrays of nanoparticles designed from fundamental structures of algebraic number theory.

I n recent years, the engineering of self-similar structures 1 in photonics and nano-optics technologies 2-7 enabled the manipulation of light states beyond periodic 8 or disordered systems 9,10 , adding novel functionalities to complex optical media [11][12][13][14] with applications to nano-devices and metamaterials [15][16][17] . The concept of multifractality (MF), which describes intertwined sets of fractals [18][19][20] , was first introduced in order to analyze multi-scale energy dissipation in turbulent fluids 20 and broadened our understanding of complex structures that appear in various fields of science and engineering 21 . Specifically, multifractal concepts provided a number of significant insights into signal analysis 22 , finance [23][24][25] , network traffic [26][27][28] , photonics 12,[29][30][31][32][33][34][35] and critical phenomena [36][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 conductors 36,40 , superconductors 38 , as well as atomic matter waves 43 . The multifractality of classical waves has also been observed in the propagation of surface acoustic waves on quasi-periodically corrugated structures 35 and in ultrasound waves through random scattering media close to the Anderson localization threshold 44 .
Considering the fundamental analogy between the behavior of electronic and optical waves 10 , 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 components 12,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 theory [46][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 characteristics 49 . 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 TiO 2 nanopillars deposited atop a transparent SiO 2 substrate and arranged according to the prime elements of the Eisenstein and Gaussian integers 47 , as well as two-dimensional cross sections of the irreducible elements of the Hurwitz and Lifschitz quaternions 48 . These arrays have recently been shown to support spatially complex resonances with critical behavior 49 akin to localized modes near the Anderson transition in random systems 40,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.
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 theorem 50 , 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) peaks 51 , an absolutely-continuous one characterized by a continuous and differentiable function 52 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 background 53,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 mathematically 51 . 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: where k is the in-plane component of the wavevector and r j are the vector positions of the N nanoparticles in the array. This is demonstrated in Fig. 1 [55][56][57][58] . Critical modes exhibit local fluctuations and spatial oscillations over multiple-length scales that are quantitatively described by multifractal analysis 45,59 . These peculiar scattering resonances have been demonstrated in one-dimensional aperiodic systems [55][56][57][58] and in the propagation of surface acoustic waves on quasi-periodically corrugated structures 35 . 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 arrays 60,61 , to visualize the dynamic of microtubules-associated protein 62 and to analyze the formation of colorimetric fingerprints of different aperiodic nano-patterned surfaces 63,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).  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 frequencyfrequency 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 ν and ν′, as defined below 65 : 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, r s = (x s , y s ) denotes the in-plane coordinates, and δIðr s ; νÞ ¼ Iðr s ; νÞ À ½ P n s¼1 Iðr s ; νÞ=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 ν ≈ 460 THz and ν ≈ 600 THz, respectively. The offdiagonal peaks of the correlation matrices reported in Figs. 3a, b In the second class, which includes the Hurwitz and Liftschitz arrays, the C(ν, ν′) 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 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(ν, ν) (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(ν, ν) 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 spectra 51,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 Qð ffiffiffiffiffiffi ffi À3 p Þ and Qð ffiffiffiffiffiffi ffi À1 p Þ, which are the commutative rings Z½ðÀ1 þ i ffiffi ffi 3 p Þ=2 and 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 ring [66][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 D f using the power-law scaling NðlÞ $ l ÀD f . Traditional (homogeneous) fractal structures are characterized by a global scale-invariance symmetry described by a single fractal dimension D f . On the other hand, heterogeneous fractals or multifractals are characterized by a continuous distribution f ðαÞ of local scaling exponents α such that NðlÞ $ l Àf ðαÞ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 multifractality 19 . 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 singularcontinuous 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) function 18 . This is reported in Supplementary Fig. 9 and discussed in the Supplementary Note 6.
A direct consequence of multifractality is the non-Gaussian nature of the probability density function (PDF) of the light intensity distribution near, or at, criticality 36,44 . We have evaluated the PDF of the scattered radiation from the histogram of the logarithm of the box-integrated intensities, i.e., Pðln I l Þ. We show in Figs. 4e-h the histograms produced by a box size of~0.8 × c/ν, where c is the speed of light and ν 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 approximation 21,44 : where D f = 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 α 0 ¼ 2μd=ð2μ þ σ 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 frequencies 44 (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 transitions 36,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 theory 69,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 waves 71,72 and stimulate the exploration of novel quantum phases of matter and many-body phenomena that emerge from fundamental structures of algebraic number theory. 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 SiO 2 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 TiO 2 thin film by reactive ion etching (RIE, Plasma-Therm, model 790) using Ar and SF 6 gases. Finally, the Cr mask was removed by wet etching in Transene 1020 73 .

Methods
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 reimages 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 linescan-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 squaresize 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(ν i , ν j ) of the frequency-frequency correlation matrices of Figs. 4a-d is the result of an average over 9 × 10 4 correlated values. Equation (2) can be written in the compact form 65 : Cðν; ν 0 Þ ¼ hδIðr; νÞδIðr; ν 0 Þi hIðr; νÞihIðr; ν 0 Þi : ð4Þ The normalization of the covariance hδIðr; νÞδIðr; ν 0 Þi with respect to the product of the average values hIðr; νÞihIðr; ν 0 Þi minimizes the intrinsic spectral effects related to the illumination source of the experimental apparatus. Each element of Cðν; ν 0 Þ 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 matrix 65  where μ(ν) is the average value and σ(ν) 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 algorithm 19 implemented in the routine FracLac (ver. September 2015) developed for the National Institutes of Health (NIH) distributed Image-J software package 74 . 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) plot 74 . 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 sizes 74 . All these aspects do not affect the main results of our paper, as shown by the error bars of Figs. 4a-d.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.