Localization landscape of optical waves in multifractal photonic membranes

: In this paper, we investigate the localization properties of optical waves in disordered systems with multifractal scattering potentials. In particular, we apply the localization landscape theory to the classical Helmholtz operator and, without solving the associated eigenproblem, show accurate predictions of localized eigenmodes for one-and two-dimensional multifractal structures. Finally, we design and fabricate nanoperforated photonic membranes in silicon nitride (SiN) and image directly their multifractal modes using leaky-mode spectroscopy in the visible spectral range. The measured data demonstrate optical resonances with multiscale intensity fluctuations in good qualitative agreement with numerical simulations. The proposed approach provides a convenient strategy to design multifractal photonic membranes, enabling rapid exploration of extended scattering structures with tailored disorder for enhanced light-matter interactions.


Introduction
In the last two decades, the search for Anderson light localization [1] and the study of disorderinduced phenomena for optical waves has stimulated the growing field of "disordered photonics", resulting in a wide range of applications to light generation and random lasing [2], solar energy [3], imaging and spectroscopy [4,5], nonlinear and quantum photonics [6][7][8].
Structurally complex photonic structures with non-periodic refractive index variations on the wavelength scale display a very rich physics driven by wave interference effects in the multiple scattering regime with profound analogies to the transport of electrons in disordered metallic alloys and semiconductors [9,10].In particular, various mesoscopic phenomena known for the electron transport in disordered materials, such as the weak localization of light [11], universal conductance fluctuations [12], and Anderson localization have found their counterparts in disordered optical materials as well [1,9].As originally understood by Philip Anderson in the context of metal-insulator transitions [1], the transport of quantum waves in strongly disordered media can be completely inhibited by exponentially localized eigenmodes of the Schrödinger equation [13].Moreover, even much weaker disorder can significantly modify the traditional Boltzmann transport picture due to recurrent scattering events that already occur in the so-called weak localization regime [11].Under these conditions, the diffusion constant is modified (i.e., renormalized) by interference phenomena giving rise to weakly-localized eigenmodes with a slow amplitude decay and large intensity fluctuations [14].
Recently, considerable progress has been made to understand the universal mechanisms of both weak and strong (Anderson) localization and efficiently predict localized modes in disordered structures based on the localization landscape theory [15][16][17].This general mathematical approach solves a Dirichlet problem defined by the uniform forcing of the system Hamiltonian and has been successfully utilized to accurately pinpoint the spatial locations of localized modes for any positive-definite disordered potential.Introducing the concept of an effective localization potential [17], researchers have been able to deterministically identify and characterize the eigenvalues and the subregions of a random medium that support localized eigenfunctions without having to solve computationally prohibitive eigenproblems.However, these current landscape-based methods are mostly focused on studying the localization of quantum waves that solve the Schrödinger equation and the question naturally arises whether the mathematical landscape theory can be extended to predict the general behavior of classical waves in complex nanophotonic structures.Recently, Anderson localization has been explored using the landscape theory in the context of classical acoustic waves by performing a Webster transformation that converts the classical wave equation into an effective Schrödinger equation with the same localization properties [18].These findings suggest a path to exploit the landscape analysis in the study of optical waves governed by the scalar Helmholtz operator, potentially enabling the rapid design and prototyping of extended photonic structures with arbitrary random potentials.In particular, this approach responds to the growing need to investigate the general localization behavior of complex optical potentials with tailored disorder beyond uncorrelated randomness.
In order to establish the potential of the localization landscape approach for the design of large-scale photonic structures with correlated disorder, we address here the optical resonances of multi-particle scattering systems with multifractal geometry.These types of random media are described by a continuous distribution of local scaling exponents that characterizes their distinctive fluctuation properties and have recently attracted a growing interest in various scientific fields ranging from finance [19], optical scattering [20] and even contemporary arts [21].They add novel functionalities to the manipulation of optical fields in complex media [20] beyond periodic [22] or uncorrelated disordered systems [7,9], with emerging applications to active nano-devices and metamaterials [23].
In this article, we propose a general method to calculate the landscape function of classical optical waves in multifractal structures based on the Helmholtz equation.First, we formulate and solve an eigenproblem analogous to the Hamiltonian operator of the Schrödinger equation and validate our results by considering the well-known problems of one-dimensional waveguide structures.We then apply the developed approach to both one-dimensional (1D) and twodimensional (2D) multifractal systems.Next, we discuss the generation of the investigated multifractal scattering structures and then we calculate the associated localization landscapes by solving the Dirichlet problem for the classical optical potential.Based on the obtained Helmholtz localization landscapes, we design and fabricate nanoperforated multifractal membranes in SiN and investigate experimentally their optical resonances using leaky-mode imaging spectroscopy.

Generation of multifractal scattering structures
In order to compute the landscape function and the eigenmodes of 1D multifractal potentials we considered two canonical models of multifractal behavior: (i) the binominal multiplicative cascade, which is based on multiplicative random processes [19,24] and (ii) the multifractal model of asset returns (MMAR), originally developed for the analysis of financial markets [19,25].
A random multiplicative cascade model is constructed based on a probability field obtained by first dividing an interval into two equal subintervals.Inside each subinterval, one assigns the probabilities   ∈ [0, 1] with  = 1, 2. This constitutes the first iteration of the process ( = 1).At the iteration step  = 2, each of the two previous subintervals is further divided in two smaller ones, and the probabilities associated with each sub-division are multiplied in random order (i.e., after random reshuffling) with the ones of the previous iterations.
At the iteration  = 3, one performs a similar division into subintervals and to each of them assigns the probabilities in random permutations from the previous iteration steps.This recurrent multiplicative cascade distribution defines a multifractal probability field in the limit of a large number of iterations [26].The probability value attached to a square region is the product of the   's of the square and of all its ancestors at previous generations.The distribution of cell values strictly depends on the initial choice of the probability vector  that acts as a control parameter for the resulting structures.The corresponding multiplicative random process is generally non-Gaussian [20].One-dimensional point patterns (i.e., point processes) with multifractal scaling properties are induced by the probability fields described above by distributing  particles on the specified line segment with probabilities that are proportional to the subinterval values.We achieved this goal using the Monte Carlo rejection scheme and we generated multifractal arrays with  = 1024 point scatterers positioned along a line segment of length  = 1024 m.
To generate the arrays, we chose the initial probability array  = [0.4,0.6].
The robustness of the proposed landscape method for the Helmholtz operator is supported by additionally considering multifractal point patterns generated based on the MMAR model [19,25].The MMAR is a continuous-time process  () that captures the heavy tails and long-memory volatility persistence of asserts often exhibited by realistic financial data [19].It is constructed by compounding a Brownian motion () with a random increasing function  () according to: where () is the price of a financial asset at a specific time and  (), known as the trading time function, is the cumulative distribution function (CDF) of a multifractal measure .In our example we compound with respect to the binomial multiplicative cascade introduced before.We computed the MMAR over a time series of length  = 1024, and assigned point scatterers along the time series using the Monte Carlo rejection method, obtaining a multifractal optical potential of a length  = 1024.Readers can find additional details on the MMAR methods in the specialized literature [19,25].

The Helmholtz localization landscape of one-dimensional multifractals
We now introduce the localization landscape approach for the Helmholtz operator in typical one-dimensional multifractal (1D) systems.The approach discussed here for random multifractal potentials can also be applied to deterministic scattering potentials with aperiodic order, which have been extensively investigated in the literature for the engineering of complex optical devices compatible with standard nanofabrication technology [6,27,28].
Our goal is to apply the localization landscape theory (LLT) to photonic problems governed by the scalar Helmholtz equation in a scattering dielectric medium with arbitrary geometry.We start from the Helmholtz equation in a non-homogeneous medium: where  0 is the free-space wavenumber and  () is the spatially varying permittivity of the medium.We then define the variation from the background medium in terms of the optical scattering potential   () as follows [6]: where  2  =  2 0   , and   is the permittivity of the background medium (we use   = 1).From the expression above we have: We can now formally express the Helmholtz equation as the equivalent Schrödinger equation: However, the potential   () can become negative when the permittivity of the scattering particles is larger than the one of the background (host medium).Therefore, in order to apply the LLT we implement the following transformation: were  max is the maximum value of the permittivity of the scattering medium.The shift introduced above renders the potential non-negative, which enables the application of the LLT to general optical wave problems without affecting their physical solutions.We now introduce the Helmholtz landscape equation based on the Dirichlet solution of the equation: where we defined the optical Hamiltonian  = −Δ +  ′ ().Equation 7allows us to investigate the localization landscape function () that encodes fundamental information on the localized eigenmodes [15][16][17].In order to numerically calculate the landscape function, we discretize the 1D potential function by considering the tight-binding problem of a single particle restricted to move along a discrete chain with first-neighbor hopping rate , on-site optical scattering potentials   , and we denote by   the landscape function value at site .Using the central difference formula [16], we can approximate Eq. 7 in the following form: where  is the discretization step of the coordinate -axis.To ensure that the spectrum of the possible eigenenergies is positive everywhere, the optical scattering potential term should satisfy the additional condition that   − 2/ 2 ≥ 0 [16].Note that since the potential  ′ is a bounded non-negative function, one can also define the so-called "effective potential"  = 1/.
The profiles of the effective potential  are often more regular than the optical potential and display clear structures of walls and wells that serve to identify the regions of localization of the eigenmodes.A rigorous connection has been made between the exponential decay of eigenfunctions and the effective potential  [17]: where  is a constant, and ℎ() is defined as the Agmon's distance [29] from  to a subset  of  defined as  = { ∈ Ω |  () ≤  + } for any  > 0, where  is the eigenvalue.This enables the accurate prediction of the geometrical supports of the localized eigenfunctions, which are centered at the locations of "wells" (i.e., local minima) of  and are bounded by the heights of the "walls" proportional to the eigenvalues [17].In order to demonstrate the predictive power of the Helmholtz landscape introduced above, we consider first the simple case of a three-layer asymmetric waveguide structure with a dielectric permittivity in the core region   =  max = 12.25 and with cladding permittivities  1 = 1.00 and  2 = 2.10.We use a wavelength  = 1.55 m and display the tight-binding optical potential from Eq. 4 and Eq. 8 in Figure 1(a).In addition, we considered a system of two coupled dielectric waveguides with the optical potential shown in Figure 1(b).These are canonical eigenproblems in photonics whose solutions can easily be obtained with traditional methods [30].Here, by computing the eigensolutions of the corresponding Schrödinger Eq. 5, we plot the first three optical modes of the two investigated systems in Figures 1(c,d), respectively.The dashed lines indicate the corresponding Helmholtz landscapes obtained by solving Eq. 7. In all the calculations we set the hopping rate to be  = 1, without loss of generality.The results plotted in Figure 1  We then apply the developed approach to the more challenging cases of fractal and multifractal optical potentials that take piecewise constant values equal to   and 1 based on whether a particle is present at a given point of the chain or is missing.The multifractal potentials generated by 1D multiplicative cascade method and by the MMAR method are shown in Figure 2(a,b), respectively.We selected the permittivity of the scattering medium to be   = 10.5 [31], and the incident wavelength  is chosen to satisfy  1 / = 0.4 where  1 is the average inter-particle separation along the chain.The selected  1 / is typical for the formation of bandgaps in photonic crystal structures [22].
The landscape functions () calculated by Eq. 7 are plotted in Figure 2(c,d).To demonstrate the accuracy of the localized mode predictions from the obtained landscape function, we directly calculated the eigenmodes from Eq. 5 of both structures and plotted some representative modes in Figure 3(a,b).Particular care has been taken to ensure that the calculated eigenvectors correspond to the eigenvalue of Eq. 5, determined by the free-space background medium, which are narrowly distributed around  2 0 (within a 5% dispersion value).We observed the coexistence of both localized modes with exponential decay in space and modes with large intensity fluctuations and reduced localization behavior as well.Moreover, in Figures 3(c,d) we compared the localized mode positions predicted by the peaks of the landscape function  pred and the actual positions  actual of the modes obtained numerically, selecting a wide range of modes that are spatially distributed along the entire line segment.Our data demonstrate almost perfect correlation between the predicted and the actual positions of the localized modes based on the computed landscape of the Helmholtz operator.

The Helmholtz localization landscape of two-dimensional multifractals
We now apply the Helmholtz landscape approach to more challenging 2D multifractal systems generated by the multiplicative cascade method as in reference [20].Therefore, we consider the natural generalization of Eq. 2: where the 2D Hamiltonian operator is defined as  = −Δ +  ′ (r).The associated localization landscape function (r) is obtained by solving the Dirichlet problem: where, similarly to the 1D case,  ′ (r) ≥ 0 everywhere in r = (, ).
The considered 2D scattering arrays are obtained from multiplicative cascade processes with initial probability vectors that correspond to monofractal (i.e., single-scaling fractals) and multifractal structures.In particular, we considered the case of an initial probability vector  = [1, 1, 1, 0] that produces a monofractal pattern and  = [1, 0.75, 0.5, 0.25] that results is a strongly inhomogeneous multifractal pattern [20].The point patterns corresponding to the generated probability random fields and obtained via Monte Carlo rejection are shown in Figures 4(a,b).Considering the typical dimensions of devices based on nanophotonic membranes, we restrict the multifractal structures within a spatial domain Ω of 50 × 50 m 2 , corresponding to approximately 2000 scattering points.The considered potentials are proportional to piecewise permittivity spatial distribution and obtained by dividing Ω into 128 × 128 unit sub-squares, and we considered a binary potential with value   to each occupied position representing the dielectric scatterers and otherwise with value 1 which represents the air.
The landscape functions are calculated by solving Eq. 11.We wish to setup the Laplacian operator with homogeneous Dirichlet boundary conditions on an open bounded domain Ω ∈ R 2 with boundary Ω.We discretize using finite elements with linear basis functions.Given a bounded, symmetric bilinear form1 (, ) that is coercive on  1 0 (Ω), we want to find  ∈  1 0 (Ω) such that  satisfies (, ) = 0, ∀ ∈  1 0 (Ω), where  1 0 (Ω) ⊂  2 (Ω), denotes the subspace of functions with square integrable derivatives that vanish on the boundary.This problem is known to have a unique solution  * [32].For the numerical comparisons in this work we consider domains that are the image of a square under a diffeomorphism, i.e., a smooth mapping from the reference domain  := [0, 1] 2 to the physical domain Ω.Quadrilateral finite element meshes and tensorized nodal basis function based on Legende-Gauss-Lobotto (LGL) points are used.We use isoparametric elements to approximate the geometry of Ω, i.e., on each element the geometry diffeomorphism is approximated using the same basis functions as the finite element approximation.The Jacobians for this transformation are computed at every quadrature point, and Gauss quadrature is used to numerically approximate integrals.We restrict our comparisons to uniformly refined conforming meshes and our implementation, written in Matlab, is publicly available [33].It does not support distributed memory parallelism, and is restricted to conforming meshes that can be mapped to a square (in 2D) or a cube (in 3D).While, in practice, matrix assembly for high-order discretizations is discouraged, we use sparse assembled operators in this prototype implementation.In order to obtain more accurate results, we oversampled the domain Ω where 2 × 2 elements per constant piece of the optical potential.We need to emphasize that the optical potentials just refers to the coordinates of scatterers and do not take into account the radii of each scatterer.That is to say we treat the scattering medium to be a point distribution, which is not the case for the fabricated nanohole membranes.However, such approximation does not affect the accuracy of the capability of the method to predict eigenmodes because the fabricated nanoholes have very small radii comparing to the multiscale nature of the scatterer separation.The multifractal potential term is added to the discrete Laplacian operator to obtain the discrete Helmholtz operator .Our numerical simulations were performed using the Boston University's Shared Computing Cluster (SCC) [34].
In our simulations, we selected the material permittivity to be   = 4, corresponding to the value for the silicon nitride (SiN) material used in the fabrication of our devices.We chose  1 / = 0.55 and  1 / = 0.62 to compute the landscape and the eigenmodes of monofractal and multifractal systems respectively.The selected  1 / lies in the typical range for a photonic crystal structure where bandgaps are likely to form [22,31] and optical modes with high-quality factors appear at the band-edges [31,35,36].We emphasize that the range of the considered spectral parameter  1 / used in the simulations overlaps with the one of the fabricated structures that will be investigated in the experimental section of this paper.
The landscape function (r) of both monofractal and multifractal structures are plotted in Figures 4(c,d).One can see that the landscape of the monofractal potential features intense  maxima localized within small groups of locally-symmetric clusters of adjacent scattering particles distributed across the structure with spatial distributions that closely resemble its fractal support.These are the regions in which the localization landscape predicts the existence of highly confined resonant modes localized by proximity effects occurring among closely coupled particles.This is the typical behavior observed in small fractal aggregates of dipolar particles [37].On the other hand, the landscape of the multifractal potential displays localization regions that are more broadly distributed across the support of the structure.This is consistent with the non-homogeneous multiscale nature of multifractal systems [20].In Figures 5(a-h) and 6(a-h), we display the calculated eigenmodes at the wavelengths specified in each panel.As for the 1D simulations, particular care has been taken to ensure that the eigenvectors are narrowly distributed around  2 0 (within a 5% dispersion value).Moreover, in order to facilitate comparisons with the experimental measurements in which multiple modes are simultaneously excited by an external source, we also plotted in Figures 5(i) and 6(i) a linear combination of the modes that independently contribute to the overall intensity distribution at  = 720 nm and  = 750 nm, respectively.As predicted by the localization landscape, we obtain strongly localized modes within small particle clusters for the monofractal structure, which are typical of self-similar  arrays where they have been successfully exploited to dramatically increase the cross sections of nonlinear and Raman processes [37].However, the situation is quite different in the multifractal structure, where high-intensity localized modes appear to spread across a more extended region of the geometrical support, reflecting its larger degree of spatial non-uniformity compared to the ones of monofractals.This is evident by looking at the mode superposition in Figure 6(i) in which the amplitude peaks cover a larger fraction of the structure.In the next sections, we discuss the fabrication and characterization of dielectric membrane structures designed based on the localization landscape theory and we demonstrate experimentally the formation of characteristic fractal modes in qualitative agreement with the landscape predictions.

Silicon nitride thin film growth and optical characterization
In this section we discuss the optical characterization of highly transparent silicon nitride thin films grown by reactive magnetron sputtering.Silicon nitride has become a well established material platform for integrated Si photonics where devices such as ultralow-loss waveguides [38], Wavelength ( high-Q resonators, and integrated structures for telecommunications, sensing, and metrology have been demonstrated [39].In this work, silicon nitride thin films were grown atop of silicon and fused silica substrates via reactive RF magnetron sputtering using a Denton Discovery 18 sputtering system with a base pressure of 2 × 10 −7 Torr and a substrate temperature held at 300 • C. We used a 99.99% purity 3-in silicon target and sputtered the thin films in a 2:1 argon-nitrogen environment at 2.5 mTorr deposition pressure and 200 W of RF power resulting in a deposition rate of ≈ 5 nm/min.All substrates were washed in a hot piranha solution and plasma ashed in an oxygen environment prior to deposition.We characterize the optical constants of the fabricated SiN thin films from the ultra-violet (UV) to the near infrared (NIR) via broadband variable angle spectroscopic ellipsometry (VASE) and normal incidence transmission measurements.In spectroscopic ellipsometry, a parameterized oscillator model is introduced to capture the relevant aspects of a materials optical dispersion as a function of wavelength.The oscillator model is then used to interpret the measured data and extract the optical constants of materials via accurate fitting procedures [40].VASE measurements at three angles (65 • ,70 • ,75 • ) were performed on SiN thin films grown atop of silicon substrates along with normal incidence transmission measurements of the same films grown atop of fused silica substrates.Figure 7 summarizes our work on characterizing the optical properties of the fabricated SiN thin films.Panels (a-b) show the measured and fitted ellipsometric parameters Ψ and Δ corresponding to the magnitude and phase of the ratio of the complex s-and p-polarized Fresnel reflection coefficients, respectively, at each angle along with each respective fit.The complex Fresnel reflection coefficients are related to the complex refractive index and thickness of the material through parameterized oscillator models.In particular, we employ a single parameterized Tauc-Lorentz (TL) oscillator model which has been used to extract the optical constants of amorphous materials in excellent agreement with experimental measurements [41].We find that this modeling method gives nearly ideal agreement with our measurements over the entire UV-VIS-NIR region investigated with a mean square error (MSE) of less than 10.In Figure 7(c) we report the extracted optical constants of the fabricated 350 nm thick SiN thin film as a function of wavelength.Furthermore, we performed normal incidence transmission measurements on the SiN thin films grown atop of fused silica substrates and independently test our TL model by plotting the expected transmission.We report our measurements and fit in fig 7(d) finding the TL model is in excellent agreement with measurement over the entire wavelength range investigated.The transmission spectra, along with the extracted optical constants of the fabricated material, demonstrates exceptionally high transparency with an almost constant refractive index of 2. Consistent with literature, we attribute these optical properties to the high RF power and substrate temperature during deposition [42].
In the next section we discuss the fabrication of nanohole perforated membranes with multifractal geometries.

Multifractal photonic membrane fabrication
We fabricate multifractal nanohole arrays within the 350 nm thick SiN thin film via electron beam lithography (EBL) and anisotropic reactive ion etching (RIE).A sketch of the process flow is shown in Figure 8(a).First a negative tone resist is spun atop of the grown material to achieve a thickness of ≈100 nm followed by soft bake.A thin, conductive coating is then spun atop of the resist to dissipate charging during the EBL process.Nano-hole patterns are then written onto the photoresist via EBL using a 30kV source (Zeiss Supra40VP).To develop the patterns, the sample is rinsed with DI water to remove the protective conductive coating and then submerged in a 4:1 ratio of developer and DI water for 55-60 seconds and blow dried with N 2 gas.A brief plasma ashing step in a pure oxygen environment is performed in order to remove residual resist left over from the development stage.This step is crucial in achieving reliable liftoff.A 20 nm chromium mask is then evaporated onto the sample (Step 3) and submerged in acetone for lift-off.A dry etch of the exposed nano-holes is performed using anisotropic reactive ion etching and finally a wet chromium etch is performed to reveal the final nano-hole structure.Representative scanning electron microscope (SEM) images of a multifractal sample are shown in Figure 8(b) along with dark-field optical images of all multifractal structures fabricated are shown in Figures 8(c-f).All patterns were written onto a square 100 m 2 area with each hole having the desired radius of 100 nm as confirmed by SEM imaging.In order to facilitate the comparison of the measured and simulated optical modes across the investigated spectral range, both monofractal and multifractal membranes were fabricated with  1 / = 0.62 at the operation wavelength  = 700 nm.

Experimental characterization of multifractal resonances
In order to demonstrate the excitation of polarization-resolved multifractal modes of the fabricated samples, we utilize a leaky mode imaging setup [43].Our experimental setup uses a broadband COMPACT K super-continuum laser source filtered by a monochromator (Oriel Cornerstone 260) to achieve a spectral linewidth of 2 nm.The filtered light is then TE polarized using a linear polarizer, collimated and focused onto the edge of the sample by an objective lens.The out of plane scattered radiation is collected by a high NA objective (50x) and an image of the modes is focused onto a camera (Thorlabs CS126MU) using a tube lens.We first measured the leaky modes of the fabricated monofractal structure as a function of wavelength from 450  nm-750 nm.In Figure 9 we show representative images of the spatial distributions of the optical modes of the monofractal structure excited using TE-polarized light at different wavelengths.We observe clearly distinct spatial distributions of the modes when excited at different wavelengths.Moreover, all mode patterns from these data feature spatial localization within small clusters of particles with highly fluctuating intensity profiles that are qualitatively similar to the numerical predictions from the landscape theory.This characteristic modal clustering behavior around smaller regions of space of the monofractal geometry is also manifested by the calculated modes in Figure 5 and is expected for self-similar structures in which scattering resonances are driven by coupling effects within locally-symmetric clusters distributed at multiple length scales [37,44].
As a comparison, we also investigated the optical resonances of the multifractal structure generated with the probability vector  = [1, 0.75, 0.5, 0.25] and a spectral parameter identified using the landscape theory, enabling a direct experimental comparison of the leaky-mode structures of fractals and multifractal systems at minimal computational cost.As shown in in Figure 10, the observed mode patterns of the multifractal structure display highly fluctuating spatial profiles distributed across the entire geometrical support, in qualitative agreement with the simulation results shown in Figure 6.Moreover, due to the spatially non-homogeneous nature of the investigated multifractal structure, the measured modes display no evident trend in the degree of their spatial localization when varying the wavelength, which is indicative of a broader spectrum of localized resonances compared to monofractals [20].We emphasize that, since the Helmholtz landscape method introduced here relies on scalar fields and point-like particles, it can only offer basic insights into the general nature of the localized modes.Specifically, these insights relate to the ability to pinpoint exactly the localization regions of the modes and predict global properties such as the localization/delocalization behavior of the resonant states.Therefore, it is unreasonable to expect a quantitative agreement with the experimental data at this stage.In particular, the model developed here does not capture the effects of the finite-size of the scattering particles, which may introduce local (small-scale) modifications to the overall mode structures compared to the ones measured on the fabricated samples.However, our work shows that the landscape approach is still valuable in providing predictive insights at low computational cost in large-scale scattering structures with complex potentials that are often beyond the reach of fully numerical, grid-based techniques such as the finite element method.Thus, we consider the proposed approach as a viable tool for the efficient first-order design of complex photonics structures.Future work will focus on the quantitative analysis of the spectral measures of the investigated multifractal optical systems, such as their integrated density of states [17], and the statistical distribution of intensity maxima of eigenvectors that unveil long-range multifractal correlations beyond the traditional Anderson model [45][46][47].

Conclusion
In this paper, we proposed a method to efficiently design extended photonic structures with multifractal geometries and investigate the fundamental properties of the fractal eigenmodes of 1D and 2D Helmholtz operators with multifractal scattering potentials.In particular, we explored two canonical multifractal systems generated from the multiplicative cascade and the MMAR methods.Without solving the associated eigenproblems, we calculated the landscape functions and accurately predicted the locations of the supported eigenmodes.Finally, based on the information obtained from the localization landscapes, we designed and fabricated multifractal photonic membranes in highly transparent SiN materials and directly imaged, using leaky-mode spectroscopy, their optical modes across the visible spectral range.The general predictions from the localization landscape of the Helmholtz operator were found to be in good qualitative agreement with the experimental data, establishing the localization landscape theory as a viable tool for the rapid exploration and benchmarking of scattering resonances in complex photonic structures with tailored multifractal disorder for nanophotonics and metamaterials applications.Finally, our results unveil the distinctive localization behavior of the scattering resonances supported by extended fractal and multifractal photonic membranes.

Fig. 1 .
Fig.1.Tight-binding optical potential of (a) asymmetric slab waveguide and (b) coupled slab waveguides system.(c,d) Spatial profiles of the lowest three optical modes from the solution of Eq. 5 along with the calculated landscape function  (doted lines) computed from Eq. 7.

Fig. 2 .
Fig. 2. Constructed multifractal potential by (a) the 1D multiplicative cascade method and (b)the MMAR method.Also shown in (c) and (d) are square root of landscape functions, normalized to their maxima corresponding to (a) and (b), calculated for  1 / = 0.4.First five modes of each potential are overlayed on top of the landscape function, normalized to their intensities respectively.

Fig. 3 .
Fig. 3. Correlation between the predicted and computed locations of the localized eigenmodes of 1D multifractal structures generated by (a) the multiplicative cascade method and (b) the MMAR method.The geometrical supports of the eigenmodes are predicted by the landscape function as explained in the text.

Fig. 4 .
Fig. 4. 2D geometries of the scattering structures corresponding to (a) the monofractal and (b) the multifractal potential.(c,d) Calculated landscape function  of the potentials shown in panels (a) and (b).The spectral parameters used in the calculations are  1 / = 0.55 and  1 / = 0.62, respectively.

Fig. 5 .
Fig. 5. (a-h) Eigenmodes whose spectral parameters vary in the range 0.55 <  1 / < 0.8.(i) Superposition of eigenmodes corresponding to a narrow range (within a 5% dispersion value) of eigenvalues around  2 0 .The wavelengths  at which the modes are computed are indicated in each panel.

Fig. 6 .
Fig. 6. (a-i) Eigenmodes whose spectral parameters vary in the range 0.6 <  1 / < 0.78.(i) Superposition of eigenmodes corresponding to a narrow range (within a 5% dispersion value) of eigenvalues around  2 0 .The wavelengths  at which the modes are computed are indicated in each panel.

Fig. 7 .
Fig. 7. Optical properties of fabricated 350nm thick SiN thin films.Ellipsometric parameters (a) Ψ and (b) Δ as measured (dashed) and fit (solid) at three angles (65 • ,70 • ,75 • ) of the deposited SiN thin films.(c) Refractive index (blue) and extinction coefficients (red) obtained from ellipsometric fitting parameters.(d) Transmission spectra at normal incidence measured (red) and fit (blue-dashed).Error bars of the measurement are represented by the red shading.

Fig. 8 .
Fig. 8. (a) Fabrication process flow via EBL and anisotropic RIE.(b) Representative SEM image of the fabricated multifractal pattern.The inset is a close up view of the same sample.The scale bars in the SEM image and inset are 4  and 600 nm respectively.(c-f) Dark field scattering images of monofractal and multifractal structures.Each array consists of ≈15,000 points written over a 100  2 area.The initial probability vector used to generate each pattern is reported above each image.

Fig. 9 .
Fig. 9. Representative leaky mode measurements of a monofractal nanohole array fabricated in SiN at various excitation wavelengths.The excitation wavelength is indicated by  in each panel.The corresponding spectral parameter varies in the range 0.6 <  1 / < 0.8.

Fig. 10 .
Fig. 10.Representative leaky mode measurements, at various wavelengths, of a multifractal nanohole array fabricated in SiN generated with an initial probability vector of p = [1, 0.75,0.5,0.25].The excitation wavelength is indicated by  in each panel.The corresponding spectral parameter varies in the range 0.6 <  1 / < 0.78.