Inverse design of metasurfaces with non-local interactions

Conventional metasurfaces have demonstrated efficient wavefront manipulation by using thick and high-aspect-ratio nanostructures in order to eliminate interactions between adjacent phase-shifter elements. Thinner-than-wavelength dielectric metasurfaces are highly desirable because they can facilitate fabrication and integration with both electronics and mechanically tunable platforms. Unfortunately, because their constitutive phase-shifter elements exhibit strong electromagnetic coupling between neighbors, the design requires a global optimization methodology that considers the non-local interactions. Here, we propose a global evolutionary optimization approach to inverse design non-local metasurfaces. The optimal designs are experimentally validated, demonstrating the highest efficiencies for the thinnest transmissive metalenses reported to-date for visible light. In a departure from conventional design methods based on the search of a library of pre-determined and independent meta-atoms, we take full advantage of the strong interactions among nanoresonators to improve the focusing efficiency of metalenses and demonstrate that efficiency improvements can be obtained by lowering the metasurface filling factors.


INTRODUCTION
Optical metasurfaces enable wavefront manipulation by spatially arranging sub-wavelength phase-shifter elements on a flat surface 1 . Almost all of the existing designs rely on the search of a library of pre-determined phase-shifter elements, meta-atoms, which introduces spatially localized phase shifts to the incident wavefront 2 . A critical requirement for the optimal design of metasurfaces is that every meta-atom must be electromagnetically independent from its neighbors 3 to assure local control of the light rather than over a macroscopic area of the metasurface. Consequently, in order to tightly localize light inside the metaatoms, most metasurfaces for wavefront control require highindex-contrast between the phase-shifter elements and the surrounding medium 4 , and fabrication of high-aspect ratio structures with deep sub-micron features and gaps, making mass volume production and hybrid integration extremely challenging [5][6][7] .
There is an increasing demand to reduce the metasurface thickness and meta-atom complexity while retaining high optical efficiency, which could circumvent the fabrication of high-aspectratio structures, and therefore facilitate higher-level integration of metasurfaces with electronics, MEMS 8,9 , and deformable substrates [10][11][12] . Recent attempts to demonstrate light manipulation with thinner waveguides have resulted in incomplete phase modulation and an intrinsic reduction in performance 13 . A promising approach is to use ultrathin dielectric nanoresonators, which exhibit electrical and magnetic dipolar resonances whose spectral response can be geometrically tuned [14][15][16] to provide both high transmission and full 2π phase control 17 . However, unlike tall waveguides with strong transversal confinement of light, ultrathin nanoresonators have long evanescent electromagnetic fields that reach beyond adjacent neighbors, coupling their resonances and making the conventional design of phase-gradient metasurfaces virtually impossible. Thus, for interacting meta-atoms, the optical properties of non-homogenous arrays, as required in a metalens layout, deviate significantly from the library of homogenous arrays. As a result, all of the thin resonator-based metasurfaces demonstrated to-date are either limited to slow and uniform phase changes, such as small-angle deflectors [18][19][20] , or limited to wavelengths longer than the visible, where materials with very high refractive index are available 21 .
On one hand, the limited availability of high-index materials in the visible spectrum represents a fundamental limitation to the design of high-efficiency ultrathin dielectric metasurfaces, mostly because the refractive indexes of all the low-absorption materials available for dielectric nanophotonics are not high enough to minimize interactions among resonators 22 . On the other hand, the conventional design approach is inherently limited in efficiency for extreme wave manipulation (e.g., large-angle deflection, nearfield focusing), due to insufficient discretization of continuous phase mapping 23 and local impedance mismatch between input and output wavefronts 24 , while strong non-local metasurfaces provide a promising solution 24 .
Inspired by recent developments in computational inverse design approaches to produce nanophotonics 25 devices and metasurfaces 26,27 with tailored responses, we have developed an inverse design strategy based on evolutionary optimization (EO) to demonstrate the thinnest (t ≈ λ/5) metalenses with the highest efficiencies in the visible spectrum. Instead of trying to avoid or minimize the interactions among meta-atoms, we build a general inverse design workflow that takes full advantage of interactions and leverages it to design high-performance metalenses. We employ a genetic algorithm (GA) to mimic natural selection in order to determine the optimal spatial arrangement of nanoresonators to achieve the desired optical performance, through an interface with finite-difference time-domain (FDTD) simulation. The optimal inverse designs are experimentally fabricated and shown to significantly improve the focusing efficiency when compared with conventional intuitive design by library search.

Ultrathin metasurfaces with non-local interactions
The studied metasurfaces are based on TiO 2 , a dielectric material fully compatible with semiconductor fabrication processes, with a relative high refractive index and negligible dissipation in the visible spectrum 28 , which has been widely used in waveguide-type metasurfaces 5,6 , but not in resonant metalenses. Our ultrathin metasurfaces consist of polarization-independent circular shaped TiO 2 resonators fabricated on a fused silica glass substrate and embedded in HSQ, a spin-on-glass material with refractive index comparable to glass and compatible with CMOS and MEMS fabrication processes 29 . The HSQ layer not only provides a protective coating for the nanostructures (similar to the SU-8 coating used in ref. 30 ), but also serves as an optical passivation layer for the resonant meta-atoms whose resonances are sensitive to the encasing environment.
In order to demonstrate the presence of strong interactions among resonators, we study the transmission and phase shifting properties of square and hexagonal arrays of the same circular resonators under illumination along its axis of symmetry (z axis) [14][15][16] . We use FDTD solvers (implemented in Lumerical) to build a database of transmission and phase values, by varying the resonator's geometrical parameters (disc radius r, edge-to-edge gap g, and thickness t). The thickness was optimized to ensure both high efficiency and robust nanofabrication. Figure 1 shows the transmission (a and d) and phase (b and e) dependence on the in-plane geometry for square (top) and hexagonal (bottom) arrays with thickness t = 115 nm. When the geometric parameters meet the Huygens condition, i.e., spectral matching of the electric dipole and magnetic dipole, the back-scattering of light by both dipoles in the incidence direction is canceled, leaving only forward propagation (white dashed circles in Fig. 1a, b, d, e) 17 . It is clearly shown that the geometrical location of the Huygens condition changes significantly with the arrangement of the nanoresonators. For the square array, the phase can be modulated while maintaining high transmission, following a path of constant pitch (white dashed lines, pitch p = 2r + g = 310 nm) through the Huygens condition region (white dashed circles). The phase shift dependence on disc radius along this path (Fig. 1c) is usually used as a library of meta-atoms to create conventional phase-gradient metasurfaces. However, the same library path (white dashed lines in Fig. 1d, e) provides significantly different responses for the hexagonal arrays (Fig. 1f). The electromagnetic fields distribution also changes with the local arrangement (insets of Fig. 1c, f), in spite of the fact that the radius and gaps between discs are exactly the same. This demonstrates that the ultrathin resonators are strongly coupled and influenced by their neighbors. As comparison, since waveguides are only weakly coupled 31 , their transmission and phase responses in the meta-atom library are almost identical for square and hexagonal arrays ( Supplementary Fig. 1). Therefore, the conventional intuitive design of matching the optical phase profile with phase-shifter elements from a library of pre-determined building blocks (Fig. 2a-c) is only applicable for independent high-aspect-ratio waveguides, but will be problematic for coupled ultrathin dielectric resonators, as the local arrangement will deviate from the pre-determined homogenous configurations.
Inverse design based on evolutionary optimization In this work, we demonstrate an inverse design strategy based on EO, which considers interactions among strongly coupled thin resonators through an interface to FDTD simulation ( Fig. 2d-f). We employ GA to mimic natural selection in order to determine the optimal arrangement of meta-atoms to achieve the target performance 32 . Our objective is to identify an arrangement of nanoresonators that maximizes the device performance i.e., the focusing efficiency. Figure 2d depicts the "ideal" axial electric energy density of a focusing metalens illuminated by a plane wave from the bottom that sends 100% of the light to a diffractionlimited focal spot at the pre-determined focal length. The inverse design workflow to achieve this target performance is summarized in Fig. 2e. In this workflow, we first randomly generate a population N p (~tens to hundreds 33 ) of metasurface structures. For each member in the population, the target properties (e.g., phase profile, transmission, and focusing intensity) are evaluated by solving Maxwell's equations (FDTD simulations) for the entire device structure. In a departure from conventional GA approach, our evolutionary procedure also takes advantage of clustering the populations into different classes; this classification is done based on performance of individual populations and the resultant grouping improves the diversity of the pool (see ref. 34 for details on the grouping). Our evolutionary process allows winning classes (representing better set of individuals) to expand their sizes; this translates into a more expanded search of the phase space allowing us to explore solutions near possible global minimum. Each class (i.e., group of population) explores a specific section of the structural phase space and avoids bias on solutions arising from use of individual populations in a purely conventional GA search. Next, our workflow computes the weighted sum of errors in prediction over all the populations, which is defined as its objective function (OF). As mentioned earlier, the OF here is set to maximize the focusing efficiency, which is computed from the intensity curve at the focal spot, as detailed in the methods section. The member I of the population is represented as a set given by Eq. (1), i ¼ x i =l x and r 0 i ¼ r i À r min =r max À r min are the normalized position and radius of meta-atom i after sorting them according to position, l x is the length of the device, r max and r min are the maximum and the minimum allowed values for the radius. The initial pool is randomly generated to start the search. The OF (focusing efficiency) values of the members in this initial pool are computed based on FDTD simulations. The members are then sorted by their OF values. A certain number of m best members among the total N p are chosen, which are then subjected to genetic operations. The probability of a member being selected for genetic operation is determined by its OF. The genetic operations included in our algorithm are (i) mutation, (ii) crossover, (iii) inversion, (iv) replacement by random structure. The details of the operation and their associated probabilities are provided in the methods section (see Supplementary Fig. 2 for a schematic representation). These operations result in "offspring" members containing the traits (binary segments) from their parent structures. Out of the N p members, and the offsprings, only the best N p set of metasurface structures are retained for the next generation. Such an optimization routine ensures that only good parameter sets survive after each generation. Upon repeating this workflow for sufficient generations, we sample diverse regions in the structural space before the run converges 35 , i.e., the OF values for a minimum members in the population are below a prescribed tolerance (see the methods section and the convergence curve in Supplementary Fig. 3). Typically, tens of such GA runs starting with different random population are necessary to scour the parameter space 36 ; we choose~10 best configurations from each of these GA runs, which give similar OF values. One such sample optimized metasurface configuration obtained using our inverse design workflow is depicted in Fig. 2f. Interestingly, we note that the typical EO designed layouts (e.g., Fig. 2f) consist of an aperiodic arrangement of nanoresonators with a geometry significantly different than the one expected from the Huygens conditions but resulting in a metalens with better performance. The inverse design workflow described above was used to optimize cylindrical metalenses for different numerical apertures (NAs) in the visible range, which were then cross-validated by experimental fabrication and characterization. In Fig. 3a, we compare the focusing efficiency η vs. NA of simulated and fabricated metalenses using both EO and conventional design methods. The main result of the study is evident in this figure: the lenses designed using EO (orange line) have focusing efficiencies much larger than conventional designs (blue line) for all the NAs considered. As expected, the efficiency decreases with increasing NA, because it is more challenging to fit the rapidly changing phase profile in the periphery of a lens with large NA. For small NAs, the phase spatial change is slow and the arrangement of phase-shifter elements does not deviate too much from the homogeneous library. The EO approach does not differ too much from the conventional design, although it slightly improves the lens focusing efficiency (e.g., NA = 0.17, η conventional~7 0%, while η EO~7 7%). For larger NAs, the phase changes more rapidly and the conventional designed layouts suffer severe deviation from the homogeneous library. In this case, the non-local interactions play an important role and the EO designed lenses significantly improve the focusing efficiency compared with conventional designs (e.g., for NA = 0.51, η conventional~3 0%, while η EO~6 0%).
In order to verify the focusing efficiency improvement of the EO designed metalenses, we fabricated lenses for green light (λ = 532 nm) with NA of 0.17, 0.51, 0.77, 0.87, and 0.92 (NA = 0.17 corresponding to a lens diameter of 12 μm and a focal length f of 35 μm, and the rest corresponding to a lens diameter of 12, 24, 36, 48 μm, respectively, with a constant f of 10 μm). We used a highprecision fabrication process, where the thickness was controlled by atomic layer deposition (ALD), and the in-plane geometry was patterned by e-beam lithography and reactive ion etching ( Supplementary Fig. 4). The experimental efficiencies are shown in Fig. 3a for the EO (orange circles) and conventional (blue circles) designs, which are in good agreement with the simulation results.
The fabricated metalenses were characterized by scanning electron microscopy (SEM) before deposition of the HSQ coating. Figure 3b (see also Supplementary Fig. 5) compares the SEMs of different lens layouts, which verifies the aperiodic nature of the EO designed lenses. The measured intensity distributions along the lens axial direction are also shown in Fig. 3b, which are in agreement with theoretical simulations (Supplementary Fig. 6). In conventional designs, the ideal phase profile cannot be implemented due to the non-local interactions between resonators and as a consequence, there is light going to diverging directions and satellite peaks, causing a reduction in the amount of light reaching the focal point. In the EO designed metalenses, the interactions among resonators are considered, resulting in an aperiodic layout with higher focusing efficiency since most of the light is directed toward the focal point. Both conventional and EO designed metalenses have diffraction-limited focal spots ( Supplementary  Fig. 7). The aperiodic layouts, especially for larger NAs, are composed of clustered discs with a larger range of size variation, leaving more empty spaces among adjacent clusters than conventional periodic layouts. As a result, the EO designs have filling factors consistently smaller than the filling factors of metalenses designed by conventional library search (Fig. 3c). The fact that a metasurface with a smaller filling factor yields higher efficiency is counterintuitive in conventional design methodology, where the phaseshifter elements have to be located as close as possible to each other in order to manipulate light efficiently 3 . This is because, instead of multiple independent phase-shifters, a group of closely spaced nanoresonators together with surrounding empty space work as a whole constitutive element, which could create complex collective optical modes to enhance or suppress scattering of light along specific directions. Similar features were recently observed in asymmetric cluster (dimer) nanoantenna arrays, which were used to improve the efficiency for large-angle deflection 23,37 . The proposed inverse computational approach introduces more possibilities for the design of metalenses while maintaining their thickness much smaller than the operating wavelength. The results of Fig. 3 indicate that the intrinsic electromagnetic interactions among nanoresonators can be used to globally engineer the metasurface response in order to optimize its farfield optical performance.

DISCUSSION
Although the metalenses designed by evolutionary algorithms demonstrated significant performance improvements, there are no reasons to assume that we have reached the performance limit for these ultrathin metasurfaces, especially for high NAs. As NA increases, the interactions play a more significant role and thus, there is more room to further improve the metalenses performance. To achieve higher efficiencies, other machine learning, optimization algorithms, and emerging concepts could be implemented, such as topology optimization which exploits multiple spatially overlapping optical modes 38,39 , and metagratings with tailored bianisotropic inclusions which enable unitary efficiency without the need of deeply sub-wavelength elements 40 .
The advent of big-data analytics has brought to the forefront powerful techniques such as deep-learning and artificial intelligence (AI)-guided approaches such as Monte Carlo Tree Search (MCTS) 41 , which may help identify better optimal configurations than those reported in this work. Although GA and other evolutionary approaches have limitations arising from slow convergence to solution, we note that the optimization of an entire spherical lens design that has non-local interactions is nontrivial for any AI or machine-learning technique. In the case of deep-learning algorithms, one typically needs to train on a data set comprising of several 100,000 to millions of focal intensity maps of lens derived from FDTD calculations. Hence, there is a lot of computational cost that is paid upfront. On the other hand, the evolutionary design requires 1000s of iterations to converge but is Fig. 2 Comparison of metalens design methodologies: conventional intuitive design by library search vs. inverse design by evolutionary optimization. a Ideal phase profile of a focusing lens. The inset equation corresponds to the 1D phase profile of a cylindrical lens, where x is the radial position relative to the lens center and f is the focal length. b A pre-determined meta-atom library (disc radius vs. phase shift diagram), similar to Fig. 1c. The phase profile is discretized using sub-wavelength meta-atoms with radius obtained from the pre-determined library, in order to generate a phase profile that is closest to the continuous phase profile. c A typical metalens layout designed by library search. The appropriate discs are placed in a constant-pitch lattice, in order to reproduce the ideal phase map. d Electric energy density distribution in the x-z cross-section of an ideal focusing lens. The lens is illuminated by an incident plane wave from the bottom and focuses the light to a focal spot at a distance f above the lens plane. e Evolutionary optimization workflow: genetic algorithm with an interface to FDTD simulation. f A typical aperiodic metalens layout designed by evolutionary optimization.
highly parallelizable. The various populations of the GA could be simultaneously run on different cores in a supercomputer. Thus, the computational efficiency to train would be high even if the cost is also high. Therefore, regardless of the optimization scheme used, the key challenge lies in the fact that non-local interactions in ultrathin metasurfaces have a complex design landscape. Navigating through this landscape necessitates the simultaneous consideration of a large number of nanoscale resonators, which is non-trivial for any deep learning, MCTS or evolutionary approach.
In the future, such computational inverse design strategies can pave the way towards realizing large-scale ultrathin metasurfaces with multiple functions and/or extreme performances at much lower fabrication demands than those being designed to-date. For example, our approach could be beneficial to recent metalenses designs involving meta-atom structures which make use of local interactions to demonstrate broadband achromatic metasurfaces 42,43 but are still limited to the search of a pre-determined library with high-aspect-ratio structures. Similarly, it could be useful to enhance recent inverse design computational frameworks that do not consider interactions among phase-shifter units 44 , in order to design metalenses for applications involving multi-wavelengths 45 or multi-angle illumination.
To the best of our knowledge, this is the first study that tries to perform a global structural optimization while concurrently accounting for non-local interactions between nanoscale resonators. This design strategy handles strongly coupled meta-atoms and enables higher efficiency metalenses with lower filling factor than the conventional designs by library search. This approach allowed us to demonstrate the highest efficiency for the thinnest (t ≈ λ/5) transmissive metalens for visible light. The low-aspectratio nanostructures are more compatible with standard massproduced fabrication techniques 46 and good candidates for building flat optical elements that can be easily integrated onto CMOS electronics and MEMS devices.

Numerical simulation
A commercial FDTD software (Lumerical) was used for nanoresonator and device simulations. To build a database, a uniform TiO 2 nanoresonator array For all the NAs studied, the EO designed metalenses are more efficient to focus light. b Scanning electron microscopy images of the meta-atoms arrangement for lenses with NA = 0.51 and 0.77, and the corresponding measured light intensity distribution along the axial direction (x-z plane) for both conventional and EO designs. c Comparative analysis of the metalens filling factors for various NAs. The metalenses designed using interactive meta-atoms have smaller fill factor than the ones designed using non-interacting ones. Size bar for the SEM is 2 μm, and 500 nm for the zoom-in.
was simulated by a single nanodisc with periodic boundary conditions in both x and y axes. The light transmission was quantified as the power transmitted through the nanoarray, normalized by the source power. The geometric parameters of both radius and gap were swept ( Fig. 1; Supplementary Fig. 1). For device simulation, a cylindrical lens was constructed by a layout of nanodiscs along the x axis, with periodic boundary conditions in the y axis, and incident wave polarization in the y axis. The focal intensity profile (|E| 2 ) was recorded and fitted to a Gaussian to calculate the focal spot power ( ffiffiffiffiffi ffi 2π p ac) and spot size (FWHM). The lens focusing efficiency was defined as the focal spot power divided by the power of a glass background with the same size as the lens.

Evolutionary algorithm
The optimization process begins by first randomly generating an initial pool of N p structures. The design constrains enforced during the search process are 0 x i l x À r i r min r i r max where g min is the minimum allowed value of the gap between adjacent meta-atoms. After the initial pool is generated, the OF for each member is computed using FDTD simulation. Then, m parent structures are chosen for genetic operations with probabilities based on the fitness, which is computed as where OF I , OF min and OF max are the OF of the member I, the minimum and maximum value of OF in the pool corresponding to the current generation. The probability for a given structure being selected is operations included in our algorithm are described as follows (Supplementary Fig. 2a).
Mutation. Mutation alters/perturbs the position and the radius of one or more meta-atoms in a given parent structure. The probability for a metaatom to be altered is randomly generated from a normal distribution.
Crossover. Crossover operations ( Supplementary Fig. 2b) involve two parents, say, parent A and parent B. In a 1 pt crossover operation, an imaginary slicing line x = x s is randomly generated. All the meta-atoms with x i < x s (x i > x s ) in parent A and x i > x s (x i < x s ) in parent B are combined to produce offspring A (offspring B). In a 2 pt crossover operation, 2 such slicing lines are generated and the meta-atoms in between the lines are swapped between parent A and parent B to produce offsprings. Similarly, in a 3 pt slicing operation 3 slicing lines are generated and the alternating segments of parent A and B are swapped to produce the offsprings.
Inversion. During inversion, the position of the meta-atoms is inverted with respect to the center of the structure.
Replacement by random structure. The offsprings during this operation is simply a completely random design. This operation helps in retaining the diversity in the pool.
Assigned probability for the genetic operations during our search were 0.6, 0.2, 0.1, 0.1 for mutation, crossover, inversion, replacement by random structure, respectively. This is partially inspired by the success of genetic algorithm for crystal structure predictions [47][48][49] , where the position of the basis atoms and the lattice vectors are optimized to have minimum energy. In addition, during the genetic evolution (especially at an early stage), some offspring designs diverge the incident light, which produce offsprings that converge light after an inversion operation. Different variations of the crossover included in our algorithm are performed with equal probability. The above-mentioned genetic operations can sometimes produce offsprings with overlapping meta-atoms. In such cases, the overlapping meta-atoms are deleted and a new meta-atom is introduced if there is enough space after satisfying the design constrains mentioned above. Out of the N p structures in the pool and the m offsprings, the best N p structures are promoted to the pool of the next generation and the cycle continues until the OFs of the best N p /8 members differ no more than 2%. The value of N p and m used during the search is 40 and 16.
All the cases were converged within 250 generations, although an improved design close to (within 3% or less than) the converged best is obtained much earlier than that, typically around 60-80 generations ( Supplementary Fig. 3). The conventional design by library search was included in the initial pool, which was not necessary but helped to speed up the EO process. Take the NA 0.51 lens as an example, starting with all random structures converged after 464 generations. The converged design is different but shares some similarity considering that the phase values are relative in an ideal phase profile ( Supplementary Fig. 8a), and also provides comparable performance (2% less in the maximum intensity).
As the device scales up, the EO process becomes computationally expensive due to increasing volume of the parameter space (number of resonators), and system size for the FDTD simulations, which scales as O(N 4 ), where N is the number of mesh grids. To improve the scalability, we simulate and optimize the entire metasurfaces, but address them in a sequence of increasing lens radius step by step, i.e., 6, 12, 18, 24 μm, with a constant f of 10 μm. In this way, for a given lens, the optimization of the central part is already converged, where the design is fixed, and GA is applied for the peripheral part only, which maintains a low volume of the parameter space and ensures efficient allocation of computational resources without repeating workload. When stitching a new peripheral part to the existing optimized central part, the boundary is offset inward to allow a 1 μm buffer zone (see an example in Supplementary Fig. 8b), which is sufficient to account the nonlocal interactions, while significantly reducing the computation expense. After optimization of the additional peripheral part, the buffer zone converges to a different design, demonstrating the effect of non-local interactions.

Nanofabrication
The metasurfaces were composed of TiO 2 nanoresonator, whose thickness was precisely controlled by ALD, and the in-plane geometry was patterned by e-beam lithography, which is shown schematically in Supplementary Fig. 4. Glass substrates with 115 nm TiO 2 ALD film were spin-coated with a bilayer of e-beam resist PMMA, and then 10 nm Au was sputtered as a conductive layer. The samples were patterned by a single step of e-beam lithography. After exposure, the Au conductive layer was removed by wet etching. The samples were developed in MIBK/IPA 3:1 at 4°C with ultrasonication for 1 min. A 10 nm Cr layer was deposited by e-beam evaporation and then lift off, forming hard masks for pattern transfer. After dry etching of TiO 2 and wet etching of Cr masks, metasurfaces of TiO 2 nanoresonator arrays were formed on the glass substrates, and then embedded in HSQ coating.

Optical characterization
An inverted microscope (Olympus IX73) was used to image the metasurface samples in transmission. A home-built set-up was used to introduce a green laser (wavelength 532 nm, Opto Engine) from the top of the microscope. The focal plane images were taken in order to characterize the focusing efficiency and focal spot size. Similar to the simulation, the measured focal intensity profile was fitted to a Gaussian distribution to calculate the focal spot intensity and FWHM. The lens focusing efficiency was defined as the focal spot intensity divided by the glass background intensity of the same size as the lens. Image stacks were taken while moving the stage in z-direction automatically (Prior ES10ZE Focus Controller) and then processed to obtain the x-z intensity distribution.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors upon reasonable request.

CODE AVAILABILITY
Code and workflow developed in this study are available from the authors upon reasonable request.