1. INTRODUCTION

The availability of increasingly powerful and affordable computational resources has opened up many new opportunities in science. Problems that do not have an analytical solution or require the evaluation of large datasets are particularly suitable to solve by computational models. The experimental examination of such problems is often also possible, but is commonly associated with high cost and time expenses. Besides that, experiments can suffer from artifacts or unexpected effects originating from perturbations in the experimental setup. For that reason, it is often desirable to have simulation models that can be used to support the experimental findings or predict the possible experimental outcomes.

In the field of photonics, a variety of established computational methods exist that can be applied on a wide range of different problems. Some of the most commonly used approaches are the finite element method (FEM), finite integral technique (FIT), finite difference time domain (FDTD) method, volume integral method (VIM), surface integral method (SIM), and fourier modal method (FMM) [1–8]. All of these methods are powerful tools that can be used to efficiently and accurately tackle several challenges in photonics. However, depending on the use case, each of these methods exhibits disadvantages that make them unsuitable for solving the problem. The FEM, FIT, FDTD, and FDFD are differential equation solvers that often lack to provide the desired efficiency due to their precision and broad generality. This is particularly a problem for systems that consist of large amounts of coupled sub-units that could be approximated by much simpler descriptions. While the FMM can efficiently deal with periodically arranged plasmonic and dielectric systems, its conventional formulation breaks down for arbitrarily arranged structures without periodicity. Thus, none of the previously mentioned methods is ideal for this study, since we aim to find an approach that efficiently predicts the scattering properties of arbitrarily distributed, oriented, and scaled plasmonic nanostructure ensembles.

Tailoring disorder has proven to be a promising approach to enhance the performance of advanced photonic applications such as solar cells [9–11] or spectrometers [12,13] and for creating and understanding structural colors [14–17]. In many cases, we can learn directly from nature [18]. Also, significant advances have been achieved in fabricating disordered structures [19–23]. Because of that, it is of special interest to further understand the influence of disorder on optical scattering properties and thus gain more control over it. Since metasurfaces exhibit widely tunable scattering properties, they hold great promise to enable many new advanced applications in photonics. Optical metasurfaces are flat arrangements of dielectric or metallic nanostructures with different sizes and orientations that allow tailoring the light propagation in a layer of subwavelength thickness [24]. Thus, it is possible to control the phase and polarization of light beams, with applications such as flat lenses [25], holograms [26], or classical and quantum sensors [27,28]. Even chiral multilayer metasurfaces with tailored disorder can be fabricated nowadays [29]. However, the initial fabrication of metasurfaces can be very time consuming, so that their experimental examination is not an effective tool for rapid prototyping. For the reasons mentioned previously, many of the existing simulation approaches are also not suitable for this task.

More sophisticated approaches have to be developed for that purpose. One option would be to describe disorder via perturbation theory. Muljarov et al. propose using the resonant-state expansion to calculate resonant states, also known as quasi-normal modes, of a perturbed system from the resonant states of an unperturbed system [30,31]. In that case, the resonant states of the unperturbed system are used as bases to set up an eigenvalue problem, which results in resonance frequencies and resonant field distributions of the perturbed system as eigenvalues and eigenvectors, respectively. The unperturbed system can be a planar waveguide [32,33] or a spherical particle [34,35], and the perturbation can be rather large if the set of basis resonant states is large enough. For instance, in Ref. [33], the resonant states of a periodic system are obtained from those of a planar slab. Alternatively, it is possible to obtain results for a metallic sphere from those of a dielectric system [34] or to switch from a spherical system to a disk [35]. While this approach is very flexible regarding geometrical modifications, it has been rarely demonstrated for modeling disordered systems consisting of several building blocks [36,37].

Alternatively, the Lemmer group in Karlsruhe has implemented a transition matrix (T-matrix) method that is optimized to be operated on graphic cards and allows for describing the interaction of thousands of spherical nanoparticles [38,39]. In that case, the analytic T-matrix for each nanoparticle is used to calculate the scattered field at each particle for a background field that consists of the incident field and the field scattered by all other spheres. While the number of nanoparticles can be rather huge, this approach has been rarely used to describe the interaction of nonspherical scatterers. In contrast, Abass et al. suggest a Green’s function method for forward and inverse modeling of quasi-periodic nanostructured surfaces [40], which is well suited for surfaces with varying height profiles.

Lalanne et al. suggest analyzing resonant states and their eigenfrequencies based on a dipolar approximation of line scatterers to investigate effects such as light localization [41]. Similarly, the Gippius group has considered the coupling of nanoantennas in periodic arrangements based on a dipolar approximation [42]. This approximation is also the basis for the well-known discrete dipole approximation [2] and used to describe surface lattice resonances [43]. We combine the various advantages of these methods into our implementation of the coupled point dipole approximation (CPDA).

This CPDA model is capable of predicting the spectrally resolved angular and spatial scattering properties of disordered plasmonic metasurfaces by using the reciprocity theorem. While we have used this model already in a recent publication to predict the visible appearance of disordered metasurfaces [16], we provide in this tutorial the theoretical background and an implementation in Python. The examined metasurfaces are perturbed with a novel disorder formalism that allows one to randomly change the location, orientation, and size of each nanostructure in a correlated or uncorrelated manner. The formalism is designed such that the disorder is progressively introduced, starting from the well-understood situation of periodically arranged nanostructures with equal dipolar properties. This enables us to slowly approach truly disordered metasurfaces and hence to analyze the effects of disorder step by step. The Cartesian and radial two-point correlation function (TPCF) as well as the structure factor (SF) are applied as metrics to quantify the disorder and relate them with the observed scattering properties of disordered metasurfaces. Together with the CPDA simulations, this provides an understanding of which scattering features are dominated by structural effects or the radiation of individual nanostructures.

The CPDA model not only provides the possibility to simulate the scattering of almost arbitrarily designed plasmonic metasurfaces, but also enables to change experimental parameters such as the NA of the objective, input and output polarization, or incidence angle of light. The latter makes the CPDA particularly suitable for computation of the highly multidimensional bidirectional reflectance distribution function (BRDF) [44]. The BRDF specifies the angular scattering of surfaces for different incidence angles, but can be extended to comprise polarization and wavelength dependence. It is regularly deployed for computer graphics or computer vision and allows one to deduce the appearance of surfaces with a dedicated color model. Vynck et al. demonstrated that it is possible to tailor the BRDF via disorder in metasurfaces. For reconstruction of the visible appearance [16,17], a modified sRGB color model can be used to compute the spatially and angle resolved appearance from the simulated hyperspectral data. This establishes the CPDA as a highly effective tool for the wavelength and polarization dependent BRDF computation of almost arbitrary plasmonic metasurfaces.

 

Fig. 1. Naturally occuring examples of surface dependent reflection. The bidirectional reflectance distribution function (BRDF) determines whether we recognize the beetles as being shiny [(a), (b) specular reflection] or matte [(c) angle-independent reflection]. Tailoring the disorder of metasurfaces can shape the BRDF to obtain the desired surface appearance and color.

Download Full Size | PDF

When looking at natural surfaces, such as those from beetles, one can see huge differences in color as well as in structural appearance. Some surfaces appear shiny with a metallic luster [see Figs. 1(a) and 1(b)], whereas others look matte [see Fig. 1(c)]. Especially, the beetles in Figs. 1(a) and 1(b) indicate that structures on the surface are extremely important for the BRDF, as they influence strongly the angle dependent reflectivity. In contrast to the beetle in (a) characterized by specular reflection, the beetle in (c) is rather dark and matte. From this simple example, one notices that tailoring the BRDF is crucial, and we ask ourselves the question whether a disordered metasurface with specific metasurface unit cells with a given three-dimensional differential scattering cross section when arranged in the appropriate fashion can tailor any desired appearance, as one would expect from rendering and ray tracing programs. Figure 2 illustrates the key question of our ansatz: can one tailor the correlation function (CF) and form factor (emission cone) of a metasurface such that any desired BRDF can be obtained?

 

Fig. 2. Concept describing the key question. (c) How can one tailor the correlation function and form factor of metasurfaces, using the individual unit cells, to obtain any desired BRDF and appearance (color, i.e., spectral reflectance amplitude, angular and polarization distribution of reflectance), as (a), (b) in the shell of beetles? The key question we study with our CPDA is what role the type and the correlation of the disorder plays for the appearance of surfaces. The images in (b) and (c) are for illustration only and do not represent the actual beetle shell.

Download Full Size | PDF

The tutorial is organized as follows: in Section 2, we introduce the basic theory related to the CPDA, the reconstruction of near and far fields via the reciprocity theorem, and disorder metrics. Furthermore, we recapitulate how to image the Fourier and image planes, and we discuss concepts for lensless imaging via Fourier transform of the far-field image and phase-retrieval algorithms. Finally, we give some introduction to color science, which is necessary for constructing the visible appearance of metasurfaces. In Section 3, we explain how our Python software package is organized. Then, we provide simulation results for the reflection of ordered and disordered metasurfaces at normal and oblique incidences in Section 4. This includes a demonstration of beam switching [45] and tailored angular scattering via the phase-retrieval algorithms described in Section 2 by our software package. Finally, we validate our model by comparison with analytical results and full-wave simulations as well as experimental results. The tutorial concludes with a summary and outlook.

2. THEORETICAL BACKGROUND

This section provides the theoretical framework behind the computational approach presented in this tutorial. All important concepts necessary to understand the content of this tutorial are summarized in three sections. In the first section, the electromagnetic interaction model describing the properties of plasmonic metasurfaces is derived. The second section introduces different types of structural disorder and provides metrics that characterize disorder. Finally, the color science behind the transition from spectral intensities to perceived colors is explained and discussed.

A. Dipolar Approximation of Plasmonic Metasurfaces

The key assumption of the introduced computational approach is that the nanostructures that constitute any plasmonic metasurface can be approximated as an ensemble of point dipoles. To understand the implications of this approximation, it is crucial to derive the most important laws starting from the fundamental Maxwell equations. This way, one can keep track of all approximations introduced and get a feeling for the strengths and limits of this approach.

The frequency domain representation of Maxwell equations is given in SI units as

Here, $\varrho$ and ${\textbf{j}}$ indicate the free charge and current densities. Furthermore, the electric displacement ${\textbf{D}}$ and magnetic induction ${\textbf{B}}$ are linked to the electric and magnetic fields ${\textbf{E}}$ and ${\textbf{H}}$, respectively, via

Assuming a spatially constant permeability $\mu$, the wave equation can be derived using Maxwell Eqs. (1c) and (1d) and Eqs. (2) and (3):

The wave equation describes an electromagnetic wave generated by a source ${\textbf{j}}$ that propagates in a medium with relative permittivity $\varepsilon$ and permeability $\mu$. Henceforth, we restrict us to spatially constant $\varepsilon$, so that we can assign the medium wavenumber $k = \sqrt {\varepsilon\mu} {k_0}$, where ${k_0} = \omega /{c_0}$ is the vacuum wavenumber.

The scalar solution of the wave equation can be found in the absence of free charges ($\nabla \cdot {\textbf{D}} = 0$) by computing the Green’s function, which is the solution for an excitation by a point source,

By using the Green’s dyadic we can also find the solution for the multidimensional wave equation driven by a point source:

It can be shown that the Green’s dyadic can be expressed in terms of the Green’s function [46]

Here, $\otimes$ denotes the outer vector product. As indicated in the second line of the equation above, the Green’s dyadic can be separated into contributions that dominate in different spatial regimes:

Although there is obviously no discrete transition between the different regimes, these contributions are commonly referred to as near field, intermediate field, and far field. The names originate from the fact that the fields dominate in different distances to the point source, depending on the ratio $1/kr$.

Knowing the Green’s dyadic, we can now find the particular solution of the wave equation in Eq. (4) by convoluting the Green’s dyadic with the driving source:

The last step in the derivation above involves the assumption that $\hat G({{\textbf{r}} - {{\textbf{r}}^\prime},\omega})$ is independent of ${{\textbf{r}}^\prime}$ for distances ${\textbf{r}}$ much larger than ${{\textbf{r}}^\prime}$, and consequently, the integral can be evaluated independently. The remaining integral can be solved by expressing the current density in terms of the dipole moment via [46]

After inserting this expression into Eq. (10), we find that the electric field radiated by a point source ${\textbf{j}}({\textbf{r}},\omega) = {\textbf{J}}(\omega)\delta ({\textbf{r}})$ is given by

By replacing the angular frequency $\omega$ with the wavenumber $k$, we can define the electric field radiated by a point dipole in terms of the wavenumber:

1. Dipolar Coupling Model

So far, we have considered only the situation of a single dipole. For a higher number of dipoles, coupling effects start to play a role and thus need to be considered in the dipole model. To describe these coupling processes mathematically, we need to recapitulate that the polarizability tensor $\hat \alpha$ is a dipole dependent property that provides the induced dipole moment ${\textbf{p}}$ for a given electric field ${{\textbf{E}}_{{\rm{tot}}}}$ at the position of that dipole. In the following derivations, the index $m = 1,2,...,N$ labels each dipole situated at position ${{\textbf{r}}_m}$ and characterized by the polarizability tensor ${\hat \alpha _m}$. This means that we can express each dipole moment induced by the total electric field ${{\textbf{E}}_{{\rm{tot}}}}$ as

Here, the total electric field consists of contributions from an external incident electric field ${{\textbf{E}}_{{\rm{inc}}}}({{{\textbf{r}}_m},k})$ and fields radiated by all other $N - 1$ dipoles ${{\textbf{E}}_{{\rm{d}},n}}({{{\textbf{r}}_m},k})$. Hence, the total electric field at position ${{\textbf{r}}_m}$ is

By inserting this expression into Eq. (14), one obtains an expression for the $m$th dipole moment in dependence of the incident electric field and all other dipole moments ${{\textbf{p}}_n}(k)$:

This set of coupled equations describes the coupling of the dipole ensemble, given an incident electric field ${{\textbf{E}}_{{\rm{inc}}}}$. By rearranging the equation, we can express our equation system as the matrix equation

At this point, it may not be entirely obvious why it is necessary to calculate the dipole moments to ultimately predict the far-field response of the dipole ensemble. However, the dipole moments are required for applying the reciprocity theorem, which is discussed in the following section.

2. Reciprocity Theorem

The reciprocity theorem forms the basis for various computational approaches in electromagnetics [41,47–50] and is also essential for the CPDA model presented in this study. Although the starting point of the mathematical derivation of the reciprocity theorem might seem less intuitive, it eventually returns the physical laws required for the CPDA model. More specifically, we follow here the idea of the near-to-far-field transformation [48,50] that is illustrated in Fig. 3 to derive the field emitted by an array of dipoles. Thus, let us consider two sources ${{\textbf{j}}_1}$ and ${{\textbf{j}}_2}$ that generate electromagnetic fields ${{\textbf{E}}_1}$ and ${{\textbf{E}}_2}$. By applying the divergence operator and cross product on these fields, we get

 

Fig. 3. Explanation of the reciprocity theorem. Within this work, the reciprocity theorem is applied for an illuminating source ${{\textbf{j}}_2}$ located at infinite distance to the sources ${{\textbf{j}}_1}$ that constitute the metasurface. The substrate beneath the sources ${{\textbf{j}}_1}$ is only for illustration and is not considered for the approach presented in this tutorial.

Download Full Size | PDF

The curl of the electric and magnetic fields can be replaced by using Maxwell equations Eqs. (1c) and (1d):

In summary, this means that we can write

We have defined the surface element ${\rm{d}}{\textbf{S}} = {\textbf{n}}{\rm{d}}S$ with the surface normal vector ${\textbf{n}}$ at the boundary of the volume $V.$ One should keep in mind that the reciprocity theorem holds true only for reciprocal materials ($\hat \varepsilon = {\hat \varepsilon ^T}$, $\hat{\mu}= {\hat{\mu}^T}$), which is the case for the materials considered in this tutorial. At this point, it is important to introduce some assumptions about our physical system to simplify the reciprocity theorem and ultimately obtain an analytical expression for the far-field response of our dipole ensemble. We define that the dipole ensemble is summarized as source ${{\textbf{j}}_1}$ and the light source that excites the dipole ensemble corresponds to source ${{\textbf{j}}_2}$. We further assume that the light source is far outside the integration volume, which means that the volume integral for source ${{\textbf{j}}_2}$ must equal zero in the considered volume, so that

If we consider a set of localized point sources as

At this point, one can already guess that the dipole moments derived in Eq. (17) are required to calculate fields ${{\textbf{E}}_1}$ and ${{\textbf{H}}_1}$ generated by the dipole ensemble. Nevertheless, we need to define our electromagnetic fields appropriately to solve the remaining surface integral. While we want to determine ${{\textbf{E}}_1}$ and ${{\textbf{H}}_1}$, we can select suitable fields ${{\textbf{E}}_2}$ and ${{\textbf{H}}_2}$ to simplify this task. Therefore, we define them to be plane waves with the mutually orthogonal transverse electric (TE) and transverse magnetic (TM) polarization vectors:

Here, $Z = \sqrt {\mu /\varepsilon}$ is the wave impedance, ${{\textbf{k}}_{\parallel}}$ is the projection of the wave vector to the $xy$ plane, i.e., it consists of wave vector components ${k_x}$ and ${k_y}$, and ${k_{\parallel}} = \sqrt {k_x^2 + k_y^2}$ is the corresponding length. Component ${k_z}$ is linked to the wavenumber via

Note that there is an ambiguity in selecting the sign of the square root. Here, we define ${k_z}$ to correspond to an outgoing plane wave [50]. Hence, the ${+}$ and ${-}$ signs in Eqs. (28a) to (28d) denote outgoing and incoming plane waves, respectively.

Since plane waves form a set of orthogonal basis functions and ${{\textbf{E}}_1}$ and ${{\textbf{H}}_1}$ are purely outgoing outside the dipole array, we may define them as a superposition of plane waves:

To find field amplitudes $\alpha _{{\rm{TE}}}^ +$ and $\alpha _{{\rm{TM}}}^ +$, we need to evaluate the integral on the left-hand side of Eq. (27). We consider planar surfaces parallel to the $xy$ plane on the top and bottom and an infinite extension to the left and right. Let us consider an incident plane wave ${{\textbf{H}}_2}$ from the top with an in-plane component ${-}{{\textbf{k}^\prime_{\parallel}}}$. We then integrate the planar surfaces at the top (incident polarization is ${{\textbf{H}}^ -}$) and bottom (incident polarization is ${{\textbf{H}}^ +}$) separately for the first term,

Here, $S$ comes from the orthonormality of plane waves and equals ${(2\pi)^2}$ for the infinitely extended surfaces. In numerical calculations, however, we truncate our computational space and assume periodic boundary conditions, so that $S$ equals the area of the surface.

This means that the complex field amplitudes depend on all dipole moments ${{\textbf{p}}_n}$ and the incident electric field ${{\textbf{E}}_2}$, which can be any arbitrary plane wave. Since the dipole moments ${{\textbf{p}}_n}$ can be computed by solving the set of equations defined in Eq. (17), we can now calculate the Fourier transform of electric field ${{\textbf{E}}_1}$, radiated into the far-field by the ensemble of point dipoles.

B. Disorder in Plasmonic Metasurfaces

With the dipole model discussed in Section 2.A, it is possible to compute the reflectance and transmittance of any arrangement of dipoles. However, to quantify the effects of disorder on plasmonic metasurfaces, it is important to introduce a disorder model that allows one to generate unbiased, randomized but still reproducible structure distributions. It is important to mention that there are many different ways to implement such a disorder model. Here, we introduce the disorder by perturbing a periodic grid of plasmonic structures. This way, it is possible to unambiguously define the non-disordered case and gradually increase the influence of the different disorder types. Figure 4 explains how positional disorder is introduced for a symmetric $3 \times 3$ grid of disks with periodicities ${p_x}$ and ${p_y}$. From now on, we refer to the number of structures as ${n_x}$ and ${n_y}$ and to the periodicities as ${p_x}$ and ${p_y}$, where $x$ and $y$ specify the corresponding spatial directions. The periodic grid is perturbed by assigning random shifts to each of the structures. These random shifts are determined by a probability density function (PDF) $P({\delta {\textbf{r}} | {s_{\!\rm{d}}}{\textbf{p}}})$, which defines the probability that a structure experiences a shift $\delta {\textbf{r}} = (\delta x, \delta y{)^T}$ in dependence of the periodicities ${\textbf{p}} = {({{p_x}, {p_y}})^T}$ and the disorder strength ${s_{\!\rm{d}}} \in {\mathbb{R}_{\geq 0}}$. For convenience, the PDF is abbreviated henceforth as $P({\delta {\textbf{r}}})$.

 

Fig. 4. Introduction of disorder and correlation. The four steps shown in this figure explain how positional disorder is introduced. The situation is similar for rotational and dimensional disorder. (1) Initially, all structures are arranged as a periodic grid with periodicities ${p_x}$ or ${p_y}$ and number of structures ${n_x}$ or ${n_y}$ in $x$ or $y$ direction, respectively. Random shift vectors (red arrows) are assigned to each structure. (2) Each structure is shifted into the randomly assigned direction. For uncorrelated disorder, the process ends here. (3) Correlation is introduced by iterating through all structures and shifting the residual structures into the same direction with the shift distance being dependent on the spatial distance to the current structure. The blue disk relates to the structure of the current iteration step, and the fading blue area represents the normal distribution that weights the shift distances. (4) Final correlated structure after repeating the iteration step in (3) for each structure.

Download Full Size | PDF

We distinguish two types of disorder: correlated and uncorrelated. For the latter, the perturbation process finishes after shifting each structure in the random direction $\delta {\textbf{r}}$. To establish a correlation between individual structure perturbations, further steps are required. First, a CF $C({\Delta {\textbf{r}}, {l_{\rm{c}}}{\textbf{p}}}) \equiv C({\Delta {\textbf{r}}})$ with the corresponding correlation length ${l_{\rm{c}}} \in {\mathbb{R}_{\geq 0}}$ needs to be introduced. This function correlates the perturbations that act on a structure at position ${{\textbf{r}}_i}$ with the perturbations acting on a structure at position ${{\textbf{r}}_j} = {{\textbf{r}}_i} + \Delta {\textbf{r}}$. The correlation is introduced by selecting one structure ${s_i}$ and shifting all other structures ${s_{\!j}}({i \ne j})$ in the same direction until each structure $i = 1, \ldots , N(N = {n_x} \cdot {n_y})$ is selected. Depending on the CF $C({\Delta {\textbf{r}}})$ and the correlation length ${l_{\rm{c}}}$, the lengths of the shift vectors are weighted differently, while the shift directions stay unaffected by the correlation process.

Thus far, uncorrelated and correlated positional disorder were explained qualitatively. However, the PDF and CF have not been specified at this point. It is important to mention that the PDF needs to be normalized by definition, which is not the case for CF. Besides that, there are only fa ew restrictions in choosing distributions for $P({\delta {\textbf{r}}})$ and $C({\Delta {\textbf{r}}})$, especially since many PDFs can be generated by applying inverse transform sampling. Within this tutorial, the uniform, normal, and triangular distributions are chosen as available distributions. The corresponding multidimensional expressions are

All three distributions are scaled such that the disorder strength and correlation length are proportional to the FWHM of the distributions, so that it is

1. Disorder Types

The previous section has given a general overview on the process of introducing disorder and correlation and also specified the underlying distributions. However, all explanations refer only to positional disorder, which affects the centroid positions of each structure. Also, rotational and dimensional disorder and any superposition of all three disorder types can be utilized to perturb the structures. All three types are shown in Fig. 5.

 

Fig. 5. Disorder types. The top panels display uncorrelated disorder, where each structure is perturbed independently of the neighboring structures. Three different geometrical perturbations of a perfect periodic array are considered here. Left: to obtain positional disorder, all structures are shifted into a random direction starting from a periodic arrangement. Middle: rotational disorder is introduced by rotating the structures around their centers. Right: dimensional disorder stretches or squeezes each structure relative to the initial structure size. In general, any superposition of these disorder types is possible. The correlated disorder in the bottom panels is different from uncorrelated disorder in respect to the dependency on the neighborhood of the introduced perturbations. More specifically, in all three cases of geometrical perturbations, the correlation depends on the spatial distance between the structures. Hence, either the random shift distances, random rotation angles, or relative size deviations are dependent on the neighboring structures.

Download Full Size | PDF

As the name indicates, rotational disorder affects the in-plane orientation of each structure. The disorder strength ${s_{\!\rm{d}}} \in {\mathbb{R}_{\geq 0}}$ is now defined relative to a full rotation of $2\pi$, which means that the PDF is now $P({\delta \phi | 2\pi {s_{\!\rm{d}}}})$. The CF is still dependent on the spatial separation $\Delta {\textbf{r}}$ of two structures but affects the rotation angle $\delta \phi$ (mathematically positive rotation with $\delta \phi = 0$ being aligned along the $x$- axis) instead of the spatial shift $\delta {\textbf{r}}$, as for positional disorder.

Similarly, dimensional disorder perturbs the size of a structure. As solely rod- and disk-shaped structures are considered within this study, we define the size as either the rod length or disk diameter, respectively. In this case, the PDF randomizes the relative size deviations $\delta S$ of a structure in dependence of the disorder strength ${s_{\!\rm{d}}} \in {\mathbb{R}_{\geq 0}}$, which yields $P({\delta S | {s_{\!\rm{d}}}})$. For instance, changing a rod with length ${S_{{\rm{pre}}}} = 200 \;{\rm{nm}}$ by $\delta S = - 0.05$ yields a rod with length ${S_{{\rm{post}}}} = ({1 + \delta S}) \cdot {S_{{\rm{pre}}}} = 190 \;{\rm{nm}}$. Again, the weights deduced from the CF depend only on the spatial separation $\Delta {\textbf{r}}$ of two structures and the correlation length ${l_{\rm{c}}}$, while now affecting the size deviations $\delta S$.

In Table 1, the dependencies of the PDFs and CFs for the three disorder types are shown. For each disorder style, separate disorder strengths, correlation lengths, and distributions can be chosen, which provides a large variety of possible disorder realizations. Since the random numbers utilized for the PDFs are generated computationally as pseudorandom numbers, it is possible to recreate a randomly disordered metasurface by specifying a so-called seed value for the pseudorandom number generator in use. Consequently, in this disorder formalism, a metasurface is fully characterized by specifying the periodicities (${p_x}$, ${p_y}$) and number of structures (${n_x}$, ${n_y}$) of the unperturbed grid as well as the disorder strengths, correlation lengths, PDFs, and CFs in addition to the random seeds.

Table 1. PDFs and CFs for Positional, Rotational, and Dimensional Disorder^a

View Table | View all tables in this article

2. Disorder Metrics

In the previous sections, we have introduced a disorder model that allows one to perturb any periodic arrangement of plasmonic structures. Thus, the question arises whether it is possible to find certain metrics that quantify the effects of disorder and correlation. At this point we investigate only the effects of positional disorder, so that the quantification of the remaining two disorder types is to be done in future work.

In this section, established metrics from astrophysics and crystallography are applied on structures generated using the disorder formalism introduced before. These metrics are the TPCF and SF and allow one to characterize metasurfaces in real and reciprocal spaces [51–53].

The TPCF is commonly used in astrophysics to characterize the non-uniformity of galaxies due to clustering. There exist various estimators that compensate for errors that originate from the spatially confined nature and finite number of structures in the investigated system [51]. A comparison of these estimators has shown that the Landy–Szalay estimator performs best in the context of astrophysics [52] and is therefore utilized to calculate the TPCFs of the disordered metasurfaces in this work. The idea behind the TPCF is to count the number of structure pairs within a certain distance interval. The Landy–Szalay estimator requires three different quantities, which are data–data pair counts $DD$, data–random pair counts $DR$, and random–random pair counts $RR$. $DD$ pair counts are calculated by iterating through all the coordinates of the structure ensemble (source) and counting the number of other structures in the source that lie within a certain distance interval. To calculate $DR$ and $RR$, it is necessary to artificially generate random coordinates used to define complete randomness. If not stated differently, the random dataset is generated by uniformly distributing coordinates in the same spatial range as the coordinates of the source. In most cases, the number of random coordinates is chosen to be $10 - 100 \times$ the number of source coordinates to ensure a good level of uniformity. Similar to $DD$, $RR$ is calculated by iterating through all coordinates of the artificial dataset and counting the number of other structures in this artificial dataset within a certain distance range. The computation of $DR$ is slightly different, as it is calculated by iterating through the source, but the number of structures in the artificial dataset within a distance interval to the source structure is counted. After computing $DD$, $DR$, and $RR$, the Landy–Szalay estimator of the TPCF is given by

Here, $F = {N_R}/{N_D}$ is the fraction of total datapoints in the random and source dataset and is required to adequately normalize the estimator. As shown in Fig. 6, we further distinguish between the radial and Cartesian TPCFs, depending on whether radial distance intervals $[{r^\prime} - \delta {r^\prime},{r^\prime} + \delta {r^\prime})$ or Cartesian intervals $[{x^\prime} - \delta {x^\prime},{x^\prime} + \delta {x^\prime}) \cap [{y^\prime} - \delta {y^\prime},{y^\prime} + \delta {y^\prime})$ are used. Both radial and Cartesian TPCFs have advantages depending on the use case. One of the greatest advantages of the radial TPCF is its one-dimensionality. This provides an easily understandable measure independent of the dimensionality of the investigated system. However, especially for highly non-isotropic systems, the radial TPCF might be deceiving. In such cases, it is inevitable to use the Cartesian TPCF for characterization of the metasurface. Besides these differences, both TPCFs have in common that positive/negative values indicate positive/negative correlation, whereas values close to zero correspond to uncorrelated randomness. This makes the TPCF a powerful metric that vividly reveals spatial correlations.

 

Fig. 6. Explanation of the radial and Cartesian TPCF. (a) Exemplary structure distribution. For the radial TPCF in (b) the data–data pair count $DD$ is calculated by counting the number of structures in the interval $[{r^\prime} - \delta {r^\prime},{r^\prime} + \delta {r^\prime})$ (green area) for any structure in the source dataset and any distance ${r^\prime}$. Similarly, the data–random pair count $DR$ and random–random pair count $RR$ can be computed, which enables to calculate the radial TPCF shown (b) using the Landy–Szalay estimator defined in Eq. (40). (c) In the case of the Cartesian TPCF, the numbers of structures in the interval $[{x^\prime} - \delta {x^\prime},{x^\prime} + \delta {x^\prime}) \cap [{y^\prime} - \delta {y^\prime},{y^\prime} + \delta {y^\prime})$ (red area) are counted. The exemplary structure distribution illustrates this for one structure in the source dataset. Repeating this for all structures in the source dataset and various distances ${x^\prime}$ and ${y^\prime}$ leads to $DD$, $DR$, and $RR$, which yields the Cartesian TPCF shown in (c). The red and green intervals shown in (a) are only for illustration and do not agree with the intervals used for the computation of the TPCFs.

Download Full Size | PDF

While the TPCF is suitable to unmask spatial characteristics, it is less convenient to unveil reciprocal properties and long-range characteristics. Therefore, we need to introduce the SF. The SF is a metric commonly used to interpret the diffraction patterns occurring in crystallography. As the name indicates, it describes the contribution to a diffraction pattern that originates from the distribution of structures. The most general form of the SF is given as

The equation above states the direct way to compute the SF of an ensemble with the structures located at positions ${{\textbf{r}}_{\textbf{n}}}$. A second possibility is to calculate the SF from Fourier transforming the Cartesian TPCF [53]. This method is commonly referred to as the Fourier transform method. However, we have found that this method is computationally more expensive and often less accurate, which is why we utilize only the direct method of Eq. (42). Besides the fact that the SF is closely related to the angular scattering intensity, it also reveals the occurrence of spatial frequencies, which makes it a powerful tool for spatial frequency analysis. This gets particularly clear when we rearrange Eq. (42):

Thus, it resembles the Fourier transform

Consequently, we have now found various metrics that are suitable to characterize the spatial and reciprocal properties of disordered metasurfaces. To get a better feeling for these metrics, we discuss the metrics for several different disordered samples in the following.

Figure 7 shows three different metrics for the unperturbed case of a periodic grid of $50 \times 50$ structures. Since this case is highly non-isotropic, the radial TPCF is less suitable for characterization. Nevertheless, the strongly pronounced and spatially very confined peaks suggest that there must be some kind of order. The Cartesian TPCF strongly resembles the unperturbed structure distribution itself. Furthermore, the spatial frequencies of ${q_x} = 2\pi /{p_x}$ and ${q_y} = 2\pi /{p_y}$ can be nicely deduced from the SF.

 

Fig. 7. Disorder metrics for the completely ordered periodic array with ${s_{\!\rm{d}}} = 0\%$ and ${l_{\rm{c}}} = 0\%$. (a) Top view on the ensemble of structures with (b)–(d) being different metrics that characterize the disorder. The black lines indicate (b) periodicity, (c) unit cell, and (d) first Brillouin zone of the unperturbed grid. Each peak in (b) and (c) corresponds to another particle in the ensemble. The decay in (b) is due to the finite size of the disordered array, which results in fewer particles at large distances. (d) The structure factor emphasizes the discrete spatial frequencies of the unperturbed structures in (a). The zeroth order is removed from the SF so that features at higher frequencies are emphasized.

Download Full Size | PDF

Figure 8 shows the same structure after increasing the disorder strength to ${s_{\!\rm{d}}} = 50\%$. The radial TPCF gains significance, since introducing disorder increases the isotropy of the ensemble. The peak positions of $\xi (r)$ are at multiples of the periodicities $p = {p_x} = {p_y}$ and exhibit a linewidth of approximately $p/2$, which is in very good agreement with the disorder strength. The same holds true for the Cartesian TPCF, which exhibits broadened peaks due to the introduced positional disorder. However, the SF still looks very much the same as the unperturbed case, which indicates that the spatial frequencies of the unperturbed sample are still dominating. Nevertheless, comparing the colorbars in Figs. 7(d) and 8(d) implies that the amplitudes have indeed dropped by a factor of two.

 

Fig. 8. Same as Fig. 7, but for ${s_d} = 50\%$ and ${l_c} = 0\%$. Both peaks exhibit broadened peaks around the periodicity of the unperturbed sample. At 50% disorder strength, particularly the diagonal peaks of the structure factor in (d) have weakened in comparison to the unperturbed case in Fig. 7. The zeroth order is removed from the SF so that features at higher frequencies are emphasized.

Download Full Size | PDF

In Fig. 9, the metasurface is depicted for a disorder strength of ${s_{\!\rm{d}}} = 200\%$. At this point, most of the features originated from the periodicity of the unperturbed grid have vanished in both TPCFs. However, both TPCFs exhibit a minimum distance below which the probability to have a neighboring structure is much lower compared to a uniformly randomized sample, which is indeed a remnant of the unperturbed periodic grid. In contrast to TPCFs, the SF still shows rather distinct peaks at the spatial frequencies of the unperturbed sample. Nevertheless, the maximum amplitude has decreased by a factor of about 50 compared to the unperturbed case in Fig. 7. This facilitates the visibility of weak spatial frequencies that arise due to the introduced disorder. These quasi-random spatial frequencies still have a lower limit at about $2.5\;{\unicode{x00B5}}{{\rm{m}}^{- 1}}$ since large spatial separations of neighboring structures are unlikely.

 

Fig. 9. Same as Fig. 7, but for ${s_d} = 200\%$ and ${l_c} = 0\%$. The peaks of the radial and Cartesian TPCFs have completely vanished at 200% disorder strength. Only a minimum distance is still visible as remnant of the initial grid. Looking closely, the structure factor indicates that more spatial frequencies are emerging in the disordered ensemble. However, spatial frequencies below $\approx 2.5\;{\unicode{x00B5}}{{\rm{m}}^{- 1}}$ are unlikely since large spatial separations between two neighboring structures are restricted for the uniform PDF.

Download Full Size | PDF

In Fig. 10, the completely disordered case is depicted, which corresponds to a disorder strength of ${s_{\!\rm{d}}} \to \infty$. In this special case, the structures are uniformly distributed over the same range as in the unperturbed case in Fig. 7. This causes an equally distributed noise to be present in all three metrics. In the case of TPCFs, this means that there is no kind of correlation, which is in good agreement with expectations. The uniform noise in the SF indicates that all spatial frequencies are equally present in the structure distribution.

 

Fig. 10. Same as Fig. 7, but for ${s_d} \to \infty$ and ${l_c} = 0\%$. The radial and Cartesian TPCFs are almost completely arbitrary in the case of a truly randomized sample. In contrast to Fig. 9, there is neither a minimum nor a maximum distance between neighboring structures. However, slightly higher values can be observed in (b) for distances below ${p_x} = {p_y}$, since the average particle density is equal to Fig. 7.

Download Full Size | PDF

The effects of correlation are shown in Fig. 11 for a correlation length of ${l_{\rm{c}}} = 400\%$ and disorder strength of ${s_{\!\rm{d}}} = 50\%$. Both TPCFs reveal that the introduced correlation has a twofold effect. The short-range correlation is increased compared to the uncorrelated case in Fig. 8. However, the long-range correlations are strongly disrupted at the same time. The radial and Cartesian TPCFs reveal that the short-range interactions reach up to two to three periodicities in each direction, which is consistent with the correlation length of ${l_{\rm{c}}} = 400\%$. These short-range correlations foster the formation of small $2 \times 2$ or $3 \times 3$ sub-grids with slightly different periodicities and orientations. This finding is reinforced by the SF, which reveals the occurrence of scattered peaks around the spatial frequencies present in the unperturbed sample in Fig. 7.

 

Fig. 11. Same as Fig. 7, but for ${s_d} = 50\%$ and ${l_c} = 400\%$. Both TPCFs in (b) and (c) show that the introduced correlation increases the short-range correlation, while breaking the long-range correlation. Particularly the Cartesian TPCF unveils that the correlation affects the four closest neighbors in $x$ and $y$ directions, which is concurrent with the correlation length of ${l_{\rm{c}}} = 400\%$. The short-range correlation facilitates the creation of many small sized sub-grids that are differently oriented and thus spread the diffraction patterns in (d) around the original positions of the unperturbed grid.

Download Full Size | PDF

C. Image- and Fourier-Plane Imaging

So far, we have discussed different metrics that allow one to characterize structural information of disordered metasurfaces. However, these metrics can be utilized only if the precise distribution of the structures is known. In the case of nanostructures, it is necessary to use diffraction-limited and high NA optics, so that the structural information can be optically resolved. However, such optical setups are usually associated with high costs. This poses the question whether there are cheaper methods that allow one to achieve diffraction-limited imaging. Indeed, there exist various approaches for so called lensless imaging that originates from x-ray applications [54,55], as the materials of conventional lenses are highly absorptive in the x-ray regime. These approaches take advantage of the Fraunhofer approximation, which states that the diffraction pattern in the far field corresponds to the Fourier transform of the diffractive element. In principle, this allows one to reconstruct the diffractive element (e.g., a metasurface) by Fourier transforming the far field. Unfortunately, due to experimental limitations, it is possible to measure only the amplitudes and not the phase of the electromagnetic far field in the x-ray and optical regimes. This problem is known as the phase problem and discards the possibility to retrieve the structure of the diffractive element via the inverse Fourier transform. However, there exist various computational approaches that make weak assumptions about the structure of the diffractive element and thus provide the possibility to reconstruct the diffractive element [56,57]. Before discussing these computational approaches, we introduce a lens-based and purely optical approach that allows one to Fourier transform and inverse Fourier transform the diffractive elements without loss of the phase information. The idea behind this approach can be understood in terms of a $4f$ system and is visualized in Fig. 12. A $4f$ system consists of two lenses (${{\rm{L}}_1}$ and ${{\rm{L}}_2}$) and a sample placed in the front focal plane of ${{\rm{L}}_1}$. Lenses ${{\rm{L}}_1}$ and ${{\rm{L}}_2}$ are arranged such that the back focal plane of ${{\rm{L}}_1}$ and the front focal plane of ${{\rm{L}}_2}$ overlap. The simplest way to quantitatively describe the $4f$ system is by using the ray transfer matrix analysis, which describes light as rays that propagate with an angle of ${\theta _x}$ at position $x$. Here, we define each ray to be described by the vector

 

Fig. 12. Explanation of the image and Fourier planes. (a) Each structure of the sample placed in the front focal plane of ${{\rm{L}}_1}$ radiates light in different directions. Two rays of the light emitted by the central structure are highlighted in red, so that they can be traced easily. All rays that originate from the same point in the front focal plane must be parallel after passing ${{\rm{L}}_1}$. This requires that these parallel beams are imaged in the back focal plane of ${{\rm{L}}_2}$. Hence, the back focal plane of ${{\rm{L}}_2}$ must be conjugated to the image plane. (b) Rays emitted in the same direction but from different locations on the sample must coincide in the back focal plane of ${{\rm{L}}_1}$. This means that each point in the back focal plane of ${{\rm{L}}_1}$ corresponds to light emitted into a certain solid angle. Since the radiation angle ${\theta _i}$ ($i = [x, y]$) relates to the reciprocal coordinates $k$ via $\sin ({{\theta _i}}) = {k_i}/k$, it is possible to convert between spatial coordinates $r$ and reciprocal coordinates $k$ with a $4f$ system for any given wavelength $\lambda = 2\pi /k$.

Download Full Size | PDF

Hence, any ray that is radiated from the front focal plane of ${{\rm{L}}_1}$ is given as

Here, $m \in \mathbb{Z}$ specifies the diffraction order. Since the spatial frequency ${q_x} = 2\pi /{p_x}$ is inversely proportional to the periodicity, we can rewrite the diffraction angle ${\theta _x}$ in terms of the spatial frequency ${q_x}$ and wavenumber $k$ as

This already proves that for any given wavelength $\lambda = 2\pi /k$, spatial frequency ${q_x}$ is proportional to propagation angle ${\theta _x}$, and hence, ${\theta _x}$ and spatial coordinate $x$ must be reciprocal quantities. Consequently, the field distribution in the back focal plane of ${{\rm{L}}_1}$ is actually the Fourier transform of the field distribution in the front focal plane of ${{\rm{L}}_1}$. If we now recall that the $x$ component of wavenumber $k$ is given by the sine of the diffraction angle, we can state that

Because of this relation, the intensity $I({\textbf{q}})$ in the Fourier plane can be decomposed into contributions from the SF $S({\textbf{q}})$ and the intensity ${| {E({\textbf{q}})} |^2}$ radiated from individual structures. Piechulla et al. [19] discussed the relation between the SF and optical quantities for disordered ensembles of PMMA nanospheres in great detail. It should not go unmentioned that the intensity distribution that results from Eq. (53b) does not consider any kind of coupling between structures.

Now that we have clarified the physical relation between the image and Fourier planes, we can move on to discuss how these insights can be utilized for lensless imaging and tailored $k$-space designs, eventually.

1. Lensless Imaging

The main idea behind lensless imaging comes from the Fraunhofer approximation, which states that the structure of a diffracting sample can be reconstructed from the diffracted wave by Fourier transforming the far field. This approach is used in various fields such as x-ray- or electron-based imaging, where the usage of lenses is difficult due to material shortcomings [54,55,58]. However, lensless imaging is also of interest for optical applications, particularly for imaging macroscopically extended samples with microscopic features that can not be placed within an optical setup.

The main limitation of this approach is that the phase and amplitude information of the electromagnetic fields in the Fourier plane must be known to reconstruct the original image. Unfortunately, the bandwidth of currently available photodetectors restricts the measurement of phase information in the optical regime. Thus, only the amplitude of optical waves is experimentally available in the Fourier plane. Nevertheless, there exist various computational approaches that allow one to recover the missing phase information under certain conditions.

Among the iterative approaches, the Gerchberg and Saxton (GS), hybrid input–output (HIO), and oversampling smoothness (OSS) algorithms are arguably the most popular and established approaches [56,57,59,60]. All of them have in common that they repeatedly transform the known Fourier-plane amplitudes into the image plane and back, while imposing certain constraints on the retrieved amplitudes and phases. To work, the available data must meet certain criteria.

First of all, it is necessary that the known Fourier-plane field amplitudes are oversampled. In other words, this means that the extent (support) of the image-plane field distribution has to be smaller than the area of the image plane that is recovered. Second, it is necessary to have a good estimate of this support region in the image plane. The size of the support region can be acquired, for instance, by low resolution optics or by using an aperture of known size that artificially limits the support.

From now on, we refer to the support region as $S$ and the electric fields in the Fourier plane as ${E_F}({\textbf{q}})$ and in the image plane as ${E_I}({\textbf{r}})$. Thus, the task of phase-retrieval algorithms is to retrieve the complex-valued fields ${E_I}({\textbf{r}})$ given the support region $S$ and the oversampled field amplitudes in the Fourier plane $| {{E_F}({\textbf{q}})} |$. The HIO and OSS algorithms tackle this problem very similarly, although there is an additional step in the OSS approach that improves the convergence, robustness against noise, and consistency [56]. In both cases, the initialization is carried out by multiplying $| {{E_F}({\textbf{q}})} |$ with a random phase. From this point, the HIO algorithm repeatedly applies the following steps.

For $i = 0$, we choose $E_I^i({\textbf{r}}) = E_I^\prime ({\textbf{r}})$ to start off. Parameter $\beta \in [{0.5, 1}]$ is the step width, and $C$ indicates the set of values that violate the image-plane constraints. These constraints $C$ can be anything known a priori about the image plane. The presumably weakest constraint is the so called nonnegativity constraint, which states that ${E_I}({\textbf{r}})$ must be real and nonnegative [60]. We have extended this constraint so that ${E_I}({\textbf{r}})$ can be chosen to be larger than a specified percentile of ${E_I}({\textbf{r}})$. This has empirically shown to be sometimes beneficial for the convergence of sparse image-plane distributions. As mentioned previously, the OSS approach is very similar but has an additional step that makes the algorithm more robust [56,57].

The function $G({{\textbf{q}}, {\sigma _i}})$ in the OSS step is a Gaussian smoothing function with zero mean and decreasing variance $\sigma _i^2$ that filters out high frequencies outside the support region $S$. A strategy that describes how ${\sigma _i}$ can be chosen is discussed in much detail in Ref. [57]. For both phase-retrieval algorithms, it is recommended to perform multiple random phase initializations to increase the chance to converge against the global minimum. We utilize the root mean squared logarithmic error (RMSLE) as a metric to judge convergence quality. For our application, it has shown to be superior to the root mean squared error (RMSE), as it lowers the influence of the strongly pronounced but less important zeroth order in the Fourier plane.

Figure 13 shows the strength of the HIO and OSS algorithms for a sample with disorder parameters ${n_x} = {n_y} = 10$, ${p_x} = {p_y} = 1\;{{\unicode{x00B5}{\rm m}}}$, ${s_{\!\rm{d}}} = 100\%$, and ${l_{\rm{c}}} = 200\%$. After only about 100 iterations and 10 initializations of the initial random phase, the HIO and OSS algorithms show already very good agreement with the actual image plane. However, the OSS algorithm further converges in contrast to the HIO algorithm and exhibits excellent agreement with the ground truth after 10,000 steps. Thus, we choose the OSS algorithm for any further results in this tutorial if not stated differently.

 

Fig. 13. Phase retrieval for lensless imaging. (a) The original structure distribution is indicated by black points. The white region is the a priori known support of the structures and specifies the approximate sample extent. (b) The structure factor of the sample in (a) is the basis of the phase-retrieval algorithm and the only known information together with the sample support $S$. (c) Comparison of the retrieved image-plane amplitudes after different numbers of steps using the HIO and OSS approaches. The OSS algorithm shows better contrast after 10,000 steps and is therefore favored.

Download Full Size | PDF

D. Color Science

Color science describes the conversion process from physical properties such as spectral intensities to perceivable colors. Thus, an accurate color model is required to predict the appearance of disordered plasmonic metasurfaces.

The CIE 1931 XYZ color space is one of the most common color spaces and allows one to convert spectral intensities to so called XYZ tristumulus values and eventually RGB values. Similar to the S-, M-, and L-cones of the human eye, the CIE standard observer XYZ color matching functions $\bar x(\lambda)$, $\bar y(\lambda)$, and $\bar z(\lambda)$ are introduced. They are shown in Fig. 14(a) and weight the spectral intensities according to the sensitivity of the human eye. For reflective or transmittive surfaces with reflectance or transmittance $P(\lambda)$, the color appearance fundamentally depends on the spectral intensity distribution $I(\lambda)$ of the illuminating light source:

 

Fig. 14. Color matching functions and CIE 1931 chromaticity diagram. (a) The CIE color matching functions $\bar x(\lambda)$, $\bar y(\lambda)$, and $\bar z(\lambda)$ weight the spectral intensities that are reflected, transmitted, or radiated by the sample of interest. By converting the obtained XYZ color space to the xyY color space, it is possible to represent any perceivable color with the chromaticity coordinates $x$ and $y$. The CIE 1931 chromaticity diagram in (b) shows all colors at a fixed luminance $Y = {Y_0}$ for different coordinates $x$ and $y.$ All visible monochrome wavelengths constitute the border of the chromaticity diagram. The corners of the triangle inside the chromaticity diagram indicate the coordinates of the red, green, and blue primaries of the sRGB gamut with the equal energy white point E lying inside the triangle. Only colors within the triangle can be represented using this sRGB gamut.

Download Full Size | PDF

However, within the frame of this work, the illuminant $I(\lambda)$ is chosen such that it has equal power at each wavelength $\lambda$. This implies that suitable reference measurements have to be conducted in the experiments to ensure that the spectral influence of the illuminating light source is cancelled out.

It is often desirable to transform the CIE XYZ color space to the CIE xyY color space, which is achieved by normalizing the XYZ tristumulus values

Since coordinate $z$ is dependent only on chromaticity coordinates $x$ and $y$, it is not required to consider it hereafter. The xyY representation allows one to distinguish between chromaticity coordinates $x$ and $y$ and luminance $Y$, which expresses only the brightness of a color independent of its chromaticity. This enables to represent any color solely by specifying chromaticity coordinates $x$ and $y$, as can be seen in the CIE 1931 chromaticity diagram in Fig. 14(b).

 

Fig. 15. Architecture of the dipolar package. The dipolar package is implemented using Python and consists of the optics, disorder, and color sub-packages. These sub-packages represent the three major fields covered in this tutorial and contain all modules in which the details of Section 2 are implemented. The optics sub-package covers all models necessary to compute light–matter interactions. Algorithms required for the generation and characterization of (disordered) metasurfaces belong to the disorder sub-package, and all color science models essential for the conversion from spectral data to RGB comprise the color sub-package.

Download Full Size | PDF

The back transformation from chromaticity coordinates $x$ and $y$ to the XYZ tristumulus can be accomplished by choosing an arbitrary luminance $Y = {Y_0}$:

The previous relations are important to understand the transition from the XYZ tristumulus to RGB values, which are connected by a single matrix ${\textbf{M}}$:

However, the entries of ${\textbf{M}}$ depend on the chosen color gamut and white point, which define the pure red, green, blue, and white of the respective color space. In contrast to the white point, the red, green, and blue primaries are adapted from the widely used sRGB gamut. The reason for this is that the equal energy illuminant E is used in this work, whereas the sRGB color space is defined in terms of the so-called D65 illuminant. Figure 14(b) and Table 2 show the coordinates of the sRGB primaries and the equal energy white point E.

Table 2. Chromaticity Coordinates of the Equal Energy White Point and sRGB Primaries

View Table | View all tables in this article

Independent of the specific RGB primaries and white point, matrix ${\textbf{M}}$ is given as

Indices $r$, $g$, $b$, and $w$ indicate that the XYZ values belong to the red, green, and blue primaries or the white point, respectively. These relations allow one to compute RGB values by inverting matrix ${\textbf{M}}$ and using Eq. (57). So far, the color model has not considered that human perception of brightness is strongly nonlinear in the physical intensity of light. This means that at this point, the RGB values are linear and require a $\gamma$ correction to match human perception. In the case of sRGB color space, the transformation from linear RGB to sRGB is defined as

3. SOFTWARE DESIGN

The aim of this section is to give a short overview on the implementation of the CPDA and the associated disorder and color science models. Thus, the CPDA package will be introduced only so far that the reader has a rough idea of its architecture and knows about the meaning of the most relevant simulation parameters and their impact on the simulation results in Section 4.

The CPDA package is named dipolar and is implemented with Python because of the availability of many well-documented open-source packages that provide the possibility to share the full package with the scientific community, once all required functionalities are implemented. A discussion about interesting functionalities that remain to be implemented can be found in Section 6.

As shown in Fig. 15, the dipolar package is divided into three sub-packages that correspond to the major fields in this tutorial, which are light–matter interaction models (optics sub-package), disorder models (disorder sub-package), and color science models (color sub-package). Together with the object-oriented implementation of many modules, this ensures the modularity and scalability of the project.

 

Fig. 16. Explanation of the CPDA simulation parameters. Despite the fact that all simulations are carried out in reflection, the simulation parameters are explained in transmission for better clarity and without loss of generality. The three bold terms correspond to the classes implemented to compute the optical response of plasmonic metasurfaces. Each instance of the Dipole class has polarizability, position, and orientation as attributes. The Ensemble class wraps all Dipole instances of a sample and allows one to access the attributes of all dipoles at once. This Ensemble instance is then passed to the Simulator class together with several simulation parameters, which are the surrounding material properties $\varepsilon$ and $\mu$, vacuum wavelength of incident light ${\lambda _0}$, polar and azimuthal angles ${\theta _{\rm{i}}}$ and ${\phi _{\rm{i}}}$ of the incident plane wave, material-independent numerical aperture of the objective ${\rm{NA}} = \sin ({{\theta _{{\max}}}})$, resolution of the spatial simulation grid ${n_{{\rm{px}},x}}$ and ${n_{{\rm{px}},y}}$, scaling factors for the simulated area in image and Fourier planes ${\delta _{\rm{i}}}$ and ${\delta _{\rm{f}}}$, and height of the simulation volume $\Delta z$. From these parameters, the complex electric field amplitudes in the Fourier plane are computed, which yield the complex electric field amplitudes in the image plane using the fourier transform (FFT). The remaining gray attributes shown above are derived from the previously mentioned parameters.

Download Full Size | PDF

A. Optics Sub-package

The optics sub-package comprises the dipole.py, ensemble.py, simulator.py, hyperspectral.py, and jones.py modules. From those, the first three modules are the most crucial, as they are inevitable for computation of the optical response of plasmonic metasurfaces. Each of these three modules contains a single class that is named like the module itself. The idea behind these classes is that each instance represents a certain unit in the real physical experiment.

An instance of the Dipole class represents an individual structure of the plasmonic metasurface and contains the polarizability, position, and orientation of that structure as attributes. The wavelength dependent polarizabilities must be known beforehand and can be computed using FEM or FIT solvers such as COMSOL Multiphysics, CST Microwave Studio, Lumerical, or any other suitable software. All polarizabilities used within this tutorial are acquired via COMSOL Multiphysics.

The Ensemble class represents the entirety of metasurfaces and contains all Dipole instances contained within this metasurface. This wrapper class provides the possibility to get and set the properties of all dipoles (e.g., position or orientation) at once and is used to construct the Simulator object that computes the optical properties of the ensemble.

In principle, it is possible to instantiate the Simulator object only by passing the corresponding Ensemble instance and the simulation wavelengths as input arguments. Figure 16 shows the most relevant simulation parameters that can be specified as attributes of the Simulator instance. These parameters are polar and azimuthal angles ${\theta _{\rm{i}}}$ and ${\phi _{\rm{i}}}$ of the incident plane wave, material independent NA of the objective ${\rm{NA}} = \sin ({{\theta _{{\max}}}})$, permittivity and permeability of the surrounding material $\varepsilon$ and $\mu$, resolution of the spatial simulation grid ${n_{{\rm{px}},x}}$ and ${n_{{\rm{px}},y}}$, scaling factors for the simulated area in image and Fourier planes ${\delta _{\rm{i}}}$ and ${\delta _{\rm{f}}}$, respectively, and height of the simulation volume $\Delta z$. Within the scope of this work, the default values of these parameters are ${\theta _{\rm{i}}}={ 0^ \circ}$, ${\phi _{\rm{i}}}={ 0^ \circ}$, ${\rm{NA}} = 1$, $\varepsilon = 2.25$, $\mu= 1$, ${n_{{\rm{px}},x}} = 256$, ${n_{{\rm{px}},y}} = 256$, ${\delta _{\rm{i}}} = {\delta _{\rm{f}}} = 1$, and $\Delta z = 0$, respectively.

Depending on these simulation parameters, the complex-valued electric field amplitudes in the Fourier plane, $\alpha _{\textit{xx}}^{\rm{f}}$, $\alpha _{\textit{xy}}^{\rm{f}}$, $\alpha _{\textit{yx}}^{\rm{f}}$, and $\alpha _{\textit{yy}}^{\rm{f}}$, are computed using the CPDA model introduced in Sections 2.A.1 and 2.A.2. Here, the subscripts indicate the polarization states of incoming or outgoing fields. With these four components, it is possible to express outgoing field ${{\textbf{E}}_{{\rm{out}}}}$ in dependence of any normalized incoming field ${{\textbf{E}}_{{\rm{in}}}}$ according to the Jones calculus as

Matrices ${\hat M_{{\rm{pre}}}}$ and ${\hat M_{{\rm{post}}}}$ are Jones matrices that allow one to manipulate the polarization state of light before and after interaction with the sample and can be calculated using the jones.py module of the optics sub-package. The superscripts and arguments of the complex field amplitudes in Eq. (61) have been left out deliberately, as it equally holds true for the complex electric field amplitudes in the image plane, $\alpha _{\textit{xx}}^{\rm{i}}$, $\alpha _{\textit{xy}}^{\rm{i}}$, $\alpha _{\textit{yx}}^{\rm{i}}$ and $\alpha _{\textit{yy}}^{\rm{i}}$, obtained via FFT. Thus, the metasurface can be probed with unpolarized light or any other polarization state that can be expressed in terms of the Jones calculus. If not stated differently, the simulation results are shown for unpolarized light.

For ${n_\lambda}$ simulation wavelengths, each of the eight complex field amplitudes $\alpha _{\textit{xx}}^{\rm{f}}$, $\alpha _{\textit{xy}}^{\rm{f}}$, $\alpha _{\textit{yx}}^{\rm{f}}$, $\alpha _{\textit{yy}}^{\rm{f}}$, $\alpha _{\textit{xx}}^{\rm{i}}$, $\alpha _{\textit{xy}}^{\rm{i}}$, $\alpha _{\textit{yx}}^{\rm{i}}$, and $\alpha _{\textit{yy}}^{\rm{i}}$ are of size ${n_{{\rm{px}},x}} \times {n_{{\rm{px}},y}} \times {n_\lambda}$. For discretization of the reciprocal space, reciprocal grid coordinates ${K_x}$ and ${K_y}$ are computed via

The auxiliary hyperspectral.py module is introduced to simplify the handling of such typically large and high dimensional datasets. Briefly speaking, the hyperspectral.py module contains the Hyperspectral class, which provides methods that allow one to represent the hyperspectral field amplitudes as spatial or spectral projection for a certain polarization state of the incoming and outgoing electric fields. Furthermore, the methods of the Hyperspectral object allow one to represent the spatial projections either as RGB appearances or as physical intensities at a certain wavelength.

B. Disorder Sub-package

In the disorder sub-package, all modules comprise those that can be used for the generation and characterization of (disordered) metasurfaces. It is important to mention that any arbitrary metasurface can be simulated using the optics sub-package and that metasurfaces are not required to be generated with the disorder sub-package. Nevertheless, the disorder sub-package particularly comes in handy if unbiased and reproducible structure distributions are required for the validation of simulation approaches, disorder metrics, or similar.

The main module of this sub-package is the disorder.py module, in which the Disorder class is implemented. This class is instantiated by passing the number of structures in $x$ and $y$ directions ${n_x}$ and ${n_y}$ as well as periodicities ${p_x}$ and ${p_y}$. By calling the positional, rotational, or dimensional methods of the Disorder instance, it is possible to introduce the different disorder types in dependence of the specified disorder strength, correlation length, PDF, CF, and random seed. It is noteworthy that the order in which the methods are called yields different results in the case of correlated disorder, since the influence of correlation depends on the spatial separation of the structures. The random numbers to create positional, rotational, and dimensional disorders are generated using the MT19937 implementation of the Mersenne-Twister pseudorandom number generator providing a specific pseudorandom sequence for a given seed [61]. If no seed is provided, an unpredictable seed is pulled from the operating system.

In the distributions.py module, the different PDFs introduced in Section 2.B.1 are implemented.

Furthermore, the metrics.py module contains the functions that compute the disorder metrics from Section 2.B.2, given a certain structure distribution.

The phase-retrieval algorithms used to recover the image-plane field distributions are implemented in the retrieval.py module and can be executed in parallel to speed up the multiple random phase initializations advised to use for better convergence.

The last module in the disorder sub-package is the interpolate.py module, which is required for introducing continuous dimensional disorder. As mentioned in the beginning of this section, the polarizability of each structure must be known previous to the CPDA simulations. However, this is a problem if the size deviations $\delta S$ are continuous, since this would require an infinite number of polarizabilities. This problem can be overcome by providing a large body of polarizabilities and interpolating the polarizabilities of intermediate-sized structures. The interpolate.py module supplies the tools required for this task and therefore provides the possibility to implement continuous dimensional disorder.

C. Color Sub-package

The color sub-package is the third and last sub-package within the dipolar package. Although cie_cmf.py is strictly speaking also a module, it actually stores only the values of the XYZ color matching functions shown in Fig. 14. These CMFs are loaded into the ColorSystem class of the color.py module and used to convert spectral intensities to the RGB color system. Thus, it is necessary to instantiate the ColorSystem object by specifying the red, green, and blue primaries as well as the white point. The subsequent conversion to the RGB color system is carried out by the methods of the ColorSystem class and strictly follows the derivations in Section 2.D.

4. SIMULATION RESULTS

The previous section briefly introducted the implementation of the CPDA and the most relevant simulation parameters required to understand the results in this section. This section does not aim to present novel metasurfaces with unique properties but rather discusses the general capabilities of the simulation model, so that the reader can assess the usability of this approach for his or her precise use case. Discussing the simulation features on a general level also provides the possibility to check the consistency of the results with fundamental physical laws. All simulation results presented in this work are computed in reflection, meaning that only backscattered light is considered. Towards the end of this section, a second approach is introduced that allows one to compute the Fourier-plane intensity distribution, and the applicability of phase-retrieval algorithms for so called $k$-space design is discussed.

A. Diffraction Limit

Conventional microscopy setups cannot distinguish two adjacent structures that are closer than a certain spatial distance. This fundamental limitation is due to the finite NA of any microscope, which thus restricts to capture the full diffraction pattern of a sample. Since the diffraction pattern comprises all structural information of the sample, the finite NA causes a limited spatial resolution. This limitation was initially described by Ernst Abbe in 1873 who found that the minimum spatial resolution $\Delta x$ depends on the ${{\rm{NA}}_{\rm{C}}}$ of the condenser, ${{\rm{NA}}_{\rm{O}}}$ of the objective, and the vacuum wavelength ${\lambda _0}$ of the light:

This general formulation of the Abbe limit can be simplified for our simulation model. First of all, in the equation above, the NA is defined as ${\rm{NA}} = n \cdot \sin ({{\theta _{{\max}}}})$. However, within this work, the NA is defined independent of the refractive index $n$ of the surrounding material and thus as ${\rm{NA}} = \sin ({{\theta _{{\max}}}})$. Second, the incident light in the CPDA model is described by a plane wave, which corresponds to central illumination, and thus, ${{\rm{NA}}_{\rm{C}}} = 0$. After adapting the diffraction limit to these specific conditions, we can reformulate the diffraction limit as

Figure 17 shows the Fourier and image planes for different NAs of the objective and a fixed wavelength ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$. All simulations are carried out for a periodic arrangement of $8 \times 8$ aluminum disks with a diameter of $0.12\;{{\unicode{x00B5}{\rm m}}}$, height of $0.06\;{{\unicode{x00B5}{\rm m}}}$, and periodicity of $0.5\;{{\unicode{x00B5}{\rm m}}}$ in both directions. The simulations vividly demonstrate the implications of a limited NA. If the NA is sufficiently large to capture all diffraction orders, it is possible to distinguish each structure in the image plane. At a NA of ${\rm{NA}} = 0.5$, the spatial resolution is $\Delta x = 0.533\;{{\unicode{x00B5}{\rm m}}}$ and therefore close to the periodicity of the sample. In this case, it is barely possible to distinguish individual structures. Only due to the finite width of the (1,0), (0,1) diffraction orders are some spatial details still present in the image plane. If the NA is further reduced to ${\rm{NA}} = 0.25$, no more spatial details are present in the image plane, since only the (0,0) diffraction order is recorded, which does not hold any spatial information. Therefore, only the borders of the sample can be resolved. Indeed, the NA not only influences the resolution limit but also affects the intensity in the image plane. This is simply caused by the circumstance that the finite NA limits the amount of backscattered optical power that can be captured by the objective. Both the NA dependent scaling of the diffraction limit and of the image-plane intensities are accurately predicted by the CPDA as shown in Fig. 17.

 

Fig. 17. NA dependent diffraction limit at normal incidence. The simulations are carried out with the CPDA for ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. (a) Strongly pronounced diffraction orders are visible in the Fourier plane due to the periodic arrangement of the dipoles. At ${\rm{NA}} = 0.5$, the (1,0) and (0,1) diffraction orders are hardly visible, resulting in a lowered resolution. At ${\rm{NA}} = 0.25$, the structural details in the image plane (b) have fully vanished, since all diffraction orders besides the (0,0) order are outside of the observed Fourier plane.

Download Full Size | PDF

B. Oblique Incidence

Besides the NA of the objective, it is also possible to change the incidence angle of the plane wave that excites plasmonic structures. This is particularly important for determination of the BRDF, for which the incidence and reflection angle must be considered. The BRDF is defined as

 

Fig. 18. Image and Fourier planes for different incidence angles. The simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.25\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. According to the law of reflection, the diffraction orders move opposed to the incidence angle in (a). This shifts some of the diffraction orders into the captured Fourier plane and consequently enhances the spatial resolution, as can be seen in (b).

Download Full Size | PDF

To get an idea of how the incidence angle affects the angle of diffraction, one can extend Eq. (49) for a reflective grating with periodicities ${p_x}$, ${p_y}$ and incidence angles of ${\theta _{{\rm{i}},x}}$, ${\theta _{{\rm{i}},y}}$ [defined as ${k_x} = k\sin ({\theta _{{\rm{i}},x}})$ and ${k_y} = k\sin ({\theta _{{\rm{i}},y}})$]:

As before, ${\theta _x}$ and ${\theta _y}$ refer to diffraction angles in $x$ and $y$ directions, respectively, while $m$ indicates the diffraction order. Since the terms on the right-hand side are constant for monochromatic light and a sample with fixed periodicity, the incidence and diffraction angles must change in opposite directions.

For better presentability, the incidence angles in Fig. 18 are specified in terms of ${\theta _{{\rm{i}},x}}$ and ${\theta _{{\rm{i}},y}}$ instead of the azimuthal and polar angles ${\phi _{\rm{i}}}$ and ${\theta _{\rm{i}}}$ shown in Fig. 16. The parameters of the sample used in Fig. 18 are chosen such that no diffraction orders are present for perpendicular incidence. It is obvious that the simulated data are self-consistent in the sense that the symmetries of the sample are also present in the Fourier- and image-plane intensity distributions for each incidence angle. Furthermore, the diffraction orders change with the different incidence angles, as stated by the law of diffraction for a reflective grating in Eq. (68). Because of that, for sufficiently high incidence angles, some of the diffraction orders move into the observable Fourier plane and thus reveal some of the previously inaccessible structural information in the image plane. As one can see, the spatial resolution increases only in the direction in which the incidence angle is changed. For that reason, the main structural details become visible at $|{\theta _{{\rm{i}},x}}| = |{\theta _{{\rm{i}},y}}| ={ 30^ \circ}$, despite sub-diffraction lattice constants ${p_x}$ and ${p_y}$.

Another interesting finding is that the maximum intensity in the image plane is enhanced for oblique incidence. This indicates that increasing the incidence angle facilitates the backscattering intensity.

C. Positional Disorder

Positional disorder describes the perturbation of the center position of each individual structure in the plasmonic metasurface. All metasurfaces shown in this section are generated using the disorder formalism introduced in Section 2.B.1. This comes with the great advantage that the performance of the CPDA model can be assessed without selection bias due to the randomized nature of the presented disorder formalism. Disorder strengths ${s_{\!\rm{d}}}$ and correlation lengths ${l_{\rm{c}}}$ are chosen as ${s_{\!\rm{d}}} = {m_{\rm{s}}} \cdot 100/3\%$ and ${l_{\rm{c}}} = {m_{\rm{c}}} \cdot 200\%$ with ${m_{\rm{s}}} \in \{0, 1, 2, 3\}$ and ${m_{\rm{c}}} \in \{0, 1, 2\}$, rerspectively. Since the correlation length has no influence on ${s_{\!\rm{d}}} = 0\%$, only ${l_{\rm{c}}} = 0\%$ is considered in that case, which yields the 10 different realizations shown in Fig. 19. It is important to mention that the size of the aluminum disks shown in Fig. 19 are to scale. This means that overlapping structures in this plot would indeed overlap in a real experimental system. The CPDA does not consider the extension of individual structures per definition and hence cannot predict the effects that originate from overlapping structures. This is particularly a problem at high disorder strengths, and low correlation lengths and can lead to severely overestimated scattering. The simulations discussed in this section are carried out for ${n_x} = {n_y} = 8$ aluminum disks with a diameter of $0.12\;{{\unicode{x00B5}{\rm m}}}$, height of $0.06\;{{\unicode{x00B5}{\rm m}}}$, and periodicities of ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$. The image- and Fourier-plane intensity distributions are evaluated at the vacuum wavelength of ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$.

 

Fig. 19. Structure and image-plane plots for positional disorder. The seed of the pseudorandom number generator to create the positional disorder is zero. The associated simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The structures are plotted to scale so that overlapping structures can be identified easily. This is particularly often the case for high disorder strengths and can cause inaccuracies in the CPDA, as the spatial extent of individual structures is not considered. Since the periodicity is above the diffraction limit, most individual structures can be identified. For the non-disordered case, the intensity at the dipole sites is almost uniform but decreases towards the borders. The uniformity of the intensities decreases with increasing disorder strengths, while the maximum intensity generally decreases for most combinations of disorder strengths and correlation lengths.

Download Full Size | PDF

Figure 19 shows that most of the disordered ensembles have decreased spatial intensities in comparison to the non-disordered case. Only at ${s_{\!\rm{d}}} = 67\%$ and ${s_{\!\rm{d}}} = 100\%$ with ${l_{\rm{c}}} = 200\%$ does a single particle show enhanced emission. This indicates that the coupling with neighboring structures is beneficial in these cases. Nevertheless, the CPDA model predicts no enhanced emission or unreasonable intensities for the two overlapping structures at ${s_{\!\rm{d}}} = 100\%$ and ${l_{\rm{c}}} = 0\%$. Thus, strong emission cannot be reached by simply placing two structures in close vicinity. Besides that, the emission of the different structures gets much less homogeneous with increasing disorder in comparison to the periodic arrangement. This agrees with expectations, since the coupling is very different, depending on the structure distribution in the proximity of each structure. As described in previous sections, only structures with a spatial separation above the diffraction limit are distinguishable in the image plane. Especially at high disorder strengths and low correlation lengths, the intensity distributions of close structures overlap.

Since Eq. (53b) states that the intensity in the Fourier plane depends on the product of the single particle radiation pattern and the SF, it is expected that the reciprocal intensities shown in Fig. 20 are dominated by the SF for unpolarized light. This means that the monochromatic intensity distributions resemble the SF and thus reveal information about the spatial frequencies and correlations that occur in that system. As expected, for the periodically arranged ensemble, only the first diffraction orders are visible in Fig. 20. Already at only ${s_{\!\rm{d}}} = 33\%$ the maximum intensity significantly decreases by a factor of about three. For increasing disorder strengths and without correlation, the diagonal (1,1) diffraction orders vanish first, until all diffraction orders are absent at ${s_{\!\rm{d}}} = 100\%$. The isotropic angular scattering at such high disorder strengths indicates that no more lattice information is present in the metasurface and hence that the structure distribution approaches a truly random distribution. If correlation is introduced into the system, it is observed that the diffraction orders are more stable against perturbations and therefore persist even at very high disorder strengths. Since the diffraction orders are scattered around their original positions for strongly disordered but correlated perturbations, it can be concluded that the introduced correlation fosters the emergence of smaller subgrids with weakly varying orientations and periodicities.

 

Fig. 20. Fourier-plane plots for positional disorder. The seed of the pseudorandom number generator to create the positional disorder is zero. The simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The (0,0) diffraction order is removed to increase the visibility of higher diffraction orders comprising the structural information. For ${s_{\!\rm{d}}} = 0\%$ and ${l_{\rm{c}}} = 0\%$, the (1,0), (0,1), and (1,1) diffraction orders are visible. In the uncorrelated case, the (1,1) diffraction orders have vanished at ${s_{\!\rm{d}}} = 67\%$, whereas the (1,0) and (0,1) modes have disappeared at ${s_{\!\rm{d}}} = 100\%$, which leads to highly isotropic scattering. Especially at intermediate disorder strengths and high correlation lengths, the diffraction orders are spread around the original locations. indicating the formation of several sub-gratings with varying periodicities and orientations.

Download Full Size | PDF

D. Rotational Disorder

Rotational disorder affects the in-plane orientation of each structure within the metasurface. Mathematically speaking, the rotation is introduced by transforming the polarizability tensor $\hat \alpha$ into a rotated coordinate system. The rotation therefore affects only structures that have distinct polarizabilities ${\alpha _{\textit{xx}}}$ and ${\alpha _{\textit{yy}}}$. Therefore, the metasurfaces presented in this section are based on rods instead of the disks in Sections 4.C and 4.E. Due to the ${C_2}$ symmetry of the rods, the disorder strength is limited to 50%, meaning that the maximum allowed rotation is ${\pm}{90^ \circ}$. The correlation lengths are chosen to be between 0% and 400%, as in Section 4.C. Besides that, the simulations were carried out for aluminum rods with size $0.12 \times 0.06 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$. For the sake of consistency with the other disorder types, the polarization state of incident light is kept unpolarized. The introduced correlation is most apparent for low disorder strengths. For higher disorder strengths, it gets increasingly difficult to spot the rotational correlation by eye. Unfortunately, no metrics have been found yet that allow one to quantify rotational disorder, so that the effects originating from correlation cannot be measured for the shown samples.

The most obvious finding for the image-plane data in Fig. 21 is that the intensities are actually enhanced for increasing rotational disorder. This suggests that the coupling between the structures is destructive in the unperturbed case and that the introduced perturbations weaken the near-field coupling and hence enhance the emission. Nevertheless, the spatial intensity variations within a metasurface and also between different metasurfaces are very weak, which indicates that the coupling effects are rather negligible for unpolarized incidence.

 

Fig. 21. Structure and image-plane plots for rotational disorder. The seed of the pseudorandom number generator to create the rotational disorder is zero. The associated simulations were carried out for ${\rm NA} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum rods with size $0.12 \times 0.06 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. Because of the ${C_2}$ symmetry of the rods, the maximum disorder strength is limited to ${s_{\!\rm{d}}} = 50\%$, which corresponds to a maximum rotation of ${\pm}{90^ \circ}$. The correlated disorder can be seen best for low disorder strengths and high correlation lengths. For any disorder strength and correlation length, the intensity distribution is very homogeneous. However, the maximum intensity is enhanced in the disordered cases in comparison to the unperturbed case. This indicates that the coupling of the dipoles is destructive, if all rods are aligned with each other.

Download Full Size | PDF

This finding is similarly true for intensity distributions in the Fourier plane shown in Fig. 22. Clearly, all angular radiation patterns are dominated by the effects of diffraction that arise from the periodic arrangement of rods. This is particularly interesting, since Eq. (53b) does not hold true anymore if the individual structures have different radiation patterns due to their distinct orientations. Therefore, it is actually necessary to describe the total radiation pattern of the metasurface with Eq. (52), which considers ensembles of non-uniform scatterers. However, the radiation patterns of the disordered surfaces strongly resemble the unperturbed case, which suggests that the field distribution of each individual scatterer is rather independent of its orientation for unpolarized incident fields.

 

Fig. 22. Fourier-plane plots for rotational disorder. The seed of the pseudorandom number generator to create the rotational disorder is zero. The simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum rods with size $0.12 \times 0.06 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. Because of the periodic arrangement of the dipole positions, the angular scattering looks very similar. Only the maximum intensity decreases by a factor of up to 1.5 due to the rotation of the dipoles.

Download Full Size | PDF

This might lead to the false conclusion that rotational effects are not suitable for so called $k$-space engineering or tailored angular scattering. Indeed, in recent years, many metasurface designs have arisen that utilize sophisticated rotational arrangements to accomplish unique spectral and angular scattering properties [45,62–64]. This means that in contrast to the findings in this section, rotational disorder can indeed be used for $k$-space engineering, if the structure orientations and incident polarizations are chosen wisely. A good approach to design rotationally disordered metasurfaces can be derived from the generalized law of reflection [63,64]:

It states that reflection angles ${\theta _x}$ and ${\theta _y}$ can be manipulated not only by incidence angles ${\theta _{{\rm{i}},x}}$ and ${\theta _{{\rm{i}},y}}$, but also by in-plane phase gradients ${\rm{d}}\Phi /{\rm{d}}x$ and ${\rm{d}}\Phi /{\rm{d}}y$. Consequently, it is possible to control the reflection angles by engineering the orientation of individual structures and thereby manipulating the in-plane phase gradients. Yin et al. [45] have introduced a metasurface based on the latter principles and exhibits customizable beam switching capabilities. These metasurfaces are discussed in Section 4.F.

E. Dimensional Disorder

The effects of size variations on individual particles can be estimated from the Mie theory for small spheres with radius $R$ and dipolar resonances. It states that the scattering cross section is proportional to the squared volume ${V_{\rm{p}}}$ of that particle [65],

 

Fig. 23. Structure and image-plane plots for dimensional disorder. The seed of the pseudorandom number generator to create the dimensional disorder is zero. The associated simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The disorder strength is limited to ${s_{\!\rm{d}}} = 50\%$ of the original size, so that unreasonably large or small diameters do not occur. Especially at high disorder strengths and correlation lengths, the size deviations are significant. The intensities are particularly strong at the locations of structures with large diameters. This finding is in good agreement with Mie theory, which states that the scattering cross section is proportional to the square of the particle volume.

Download Full Size | PDF

Similar to previous disorder types, the metasurfaces in this section consist of ${n_x} = {n_y} = 8$ aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$ and periodicities ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$. The intensity distributions shown in Figs. 23 and 24 are evaluated for ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$. As stated in Section 2, only the particle diameter is affected by the dimensional disorder. The disorder strength is limited at ${s_{\!\rm{d}}} = 50\%$ to avoid unreasonably small or even negative particle diameters. Since the volume of a disk scales with the square of the diameter, it is expected that the spatial intensity scales with the fourth power of its diameter.

 

Fig. 24. Fourier-plane plots for dimensional disorder. The seed of the pseudorandom number generator to create the dimensional disorder is zero. The simulations were carried out for ${\rm{NA}} = 1$, ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. Similar to rotational disorder, the diffraction orders appear to be rather unaffected by the dimensional disorder. However, particularly at ${s_{\!\rm{d}}} = 50\%$ and ${l_{\rm{c}}} = 400\%$, the diffraction orders get blurry due to the strongly varying dipole sizes. Additionally, the angular backscattering is intensified as well, if the sample consists of many large disks with high scattering cross sections.

Download Full Size | PDF

In fact, Fig. 23 indicates that this scaling behavior is in good agreement with Eq. (70). For instance, at ${s_{\!\rm{d}}} = 50\%$ and ${l_{\rm{c}}} = 400\%$, the largest disks are approximately twice as large as in the unperturbed case, causing a $286.7/16.8 \approx 17 \times$ enhanced emission. This is very close to the theoretically expected $16 \times$ enhancement, especially regarding the distinct particle shape and the lack of coupling considerations in Eq. (70). The maximum intensity is enhanced for any disordered metasurface because the dimensional disorder always causes some disks to increase in diameter in the shown examples. Consequently, the highest spatial intensity is predominantly determined by the size of the largest disk in the ensemble.

Similar to the findings for rotational disorder, the angular scattering behavior of the dimensionally disordered metasurfaces depicted in Fig. 23 is rather dominated by the periodic arrangement of the structures. However, the amplitude variations between the different Fourier-plane intensities are significantly different between the various metasurfaces. In contrast to the image-plane intensities, the Fourier-plane intensities are enhanced only by a factor of $150.3/53.4 \approx 3$ between the unperturbed surface and disordered surface at ${s_{\!\rm{d}}} = 50\%$ and ${l_{\rm{c}}} = 400\%$. This is primarily for two reasons. First of all, the Fourier-plane intensities correspond to the differential scattering cross section

F. Beam Switching

Yin et al. [45] have introduced a metasurface that allows one to switch the deflection angle of a circularly polarized beam. The presented metasurface comprises metallic rods stepwise rotated along the $x$ axis, which causes a constant in-plane phase gradient. The rotation steps are chosen such that a full rotation from zero to $\pi$ is accomplished within one super-period $P = {n_{\rm{s}}} \cdot p$. In this notation, ${n_{\rm{s}}}$ refers to the number of structures within one super-period and $p = {p_x} = {p_y}$ to the spatial separation of the structures in $x$ and $y$ directions. The generalized law of reflection for this metasurface can thus be derived from Eq. (69) for perpendicular incidence:

Obviously, the deflection angle ${\theta _x}$ only depends only on the wavelength $\lambda = {\lambda _0}/n$ of the incident field and super-period $P.$ To investigate the influence of the total number of structures, we introduce a third parameter ${n_{\rm{p}}}$ that indicates the number of super-periods of the metasurface.

Metasurfaces are designed such that the number of structures in both directions is given by ${n_x} = {n_y} = {n_{\rm{s}}} \cdot {n_{\rm{p}}}$. This provides the possibility to examine the influence of the interstructure distance $p$, number of structures within one super-period ${n_{\rm{s}}}$, number of super-periods ${n_{\rm{p}}}$, and thus, also the total number of particles $N = {n_x} \cdot {n_y} = ({n_{\rm{s}}} \cdot {n_{\rm{p}}})^2$ independently.

The directivity

Reducing the distance $p$ between structures has essentially two effects. First of all, for constant ${n_{\rm{s}}}$, the super-period $P$ is equally decreased, which increases the angle of deflection. For the same reason, the angle of deflection stays unchanged if ${n_{\rm{s}}}$ and $p$ are changed such that the super-period $P = {n_{\rm{s}}} \cdot p$ stays unaffected.

The simulated data in Figs. 25(b) and 25(c) are based on the same type of metasurface for ${n_{\rm{s}}} = {n_{\rm{p}}} = 4$, ${p_x} = {p_y} = p = 0.3\;{{\unicode{x00B5}{\rm m}}}$, and aluminum rods with size $0.12 \times 0.06 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. As previously, the depicted Fourier-plane intensities are evaluated at ${\lambda _0} = 0.6\;{{\unicode{x00B5}{\rm m}}}$, but in contrast to the previous simulation, an additional read-out polarization is specified. Yin et al. [45] have shown that the polarization of the deflected beam is changed from right-handed circularly polarized (RCP) to left-handed circularly polarized (LCP), if the beam is transmitted. However, in Fig. 25, the deflected beam keeps its polarization state. This originates from the fact that reflection at an interface reverses the handedness of circularly polarized light. For that reason, the light reflected into the zeroth order changes the polarization state from LCP to RCP or vice versa. We also find that the sign of the deflection angle ${\theta _x}$ depends on the circular polarization state of the incident field. Consequently, both diffraction orders are present if the incident beam is unpolarized.

 

Fig. 25. Parameters and polarization dependence of the beam switching metasurface. (a) The beam switching is accomplished by metasurfaces that consist of metallic rods rotated around their center dependent on the position along the $x$ axis. Since the distance between two neighboring rods is fixed at ${p_x} = {p_y} = p$, the diffraction angle is solely determined by the number of structures ${n_{\rm{s}}}$ within one of the ${n_{\rm{p}}}$ super-periods $P = {n_{\rm{s}}} \cdot p$. All metasurfaces in this section have an equal number of structures in both directions ${n_x} = {n_y} = {n_{\rm{s}}} \cdot {n_{\rm{p}}}$. The simulations in (b) and (c) were carried out for ${\lambda _0} = 0.6\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = {n_{\rm{s}}} \cdot {n_{\rm{p}}}$, ${n_{\rm{s}}} = {n_{\rm{p}}} = 4$, ${p_x} = {p_y} = p = 0.3\;{{\unicode{x00B5}{\rm m}}}$, and aluminum rods with size $0.12 \times 0.06 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. (b) Both diffraction orders are present in the Fourier plane for unpolarized illumination and no manipulation of the polarization state at the output. (c) Incident light with circular polarization is either diffracted or stays unaffected by the metasurface depending on the read-out polarization state at the output.

Download Full Size | PDF

G. Structure Factor Approximation

The structure defined in Eq. (41) not only can be used to characterize disorder in metasurfaces, but also provides the possibility to predict Fourier-plane intensity. Although Eq. (52) allows one to compute intensity distributions for arbitrary metasurfaces, we limit our discussion to surfaces with equal structures, so that the simplified formulation in Eq. (53b) applies. Piechulla et al. [19,66] have discussed this approach in great detail for nearly hyperuniform metasurfaces that suppress long-range density fluctuations and thus promise to eliminate the zeroth order of the SF [67–69]. We refer to the approach derived from Eq. (53b) as structure factor approximation (SFA). As depicted in Fig. 26, the scattering amplitudes of individual aluminum rods are obtained from the CPDA approach for incident light polarized along the long axis of the rod. Since the SF does not depend on the shape of the individual structures, it is possible to separate the total response of the metasurface into contributions of the individual structures and the distributions of those.

 

Fig. 26. Explanation of the structure factor approximation (SFA) (top) and comparison of the SFA and CPDA (bottom). The seed of the pseudorandom number generator to create the positional disorder is zero. The simulations were carried out for ${\lambda _0} = 0.6\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$, ${s_{\!\rm{d}}} = 100\%$, ${l_{\rm{c}}} = 400\%$, and aluminum rods with size $0.18 \times 0.09 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The single particle intensity is computed for $x$-polarized light using the CPDA. Here, the structure factor (SF) does not consider the shape of individual particles and therefore returns only the spatial frequency distribution of the dipole positions. Since the sample consists of rods with equal polarizability, it is possible to compute the Fourier-plane intensity distribution by multiplying the single particle intensity with the SF via the SFA. Note that the polarization of incident light is aligned with the rod orientation. As seen in the bottom panels, the agreement of the CPDA and SFA is very good and mainly differs in the maximum predicted intensities, which is slightly lower for the CPDA. This can be explained by the coupling between the dipoles, which is not considered in the SFA.

Download Full Size | PDF

Figure 26 vividly illustrates that the envelope of the total Fourier-plane response is determined by the individual structures, while the integral subtleties arise due to the global structure distribution comprising the SF. This comes from the reciprocal relation of ${\textbf{r}} = (x, y{)^T}$ and ${\textbf{q}} = ({q_x}, {q_y}{)^T}$, which imposes an inverse relationship of the spatial extents in both spaces.

Figure 26 compares the SFA with the CPDA for the same structure distribution but different rod orientations. The polarization of incident light is always aligned with the orientation of aluminum rods. This means that the SF stays unaffected while the scattering amplitudes of the individual structures are altered corresponding to the rod orientation. Both computational approaches yield angular radiation patterns that resemble the SF with different regions of it being more or less pronounced depending on the orientation of the rods. Generally, the agreement of the SFA and CPDA is very good, with only very subtle differences. Figure 26 demonstrates that both models mainly differ in the maximum predicted intensities for which the CPDA predicts slightly lower values than the SFA. It is likely that this is due to the lack of coupling effects in the SFA, especially since Fig. 22 also indicates that the coupling induces weakened emission for aluminum disks and distances around ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$. Because of that, it is expected that the differences between the CPDA and SF are greater for metasurfaces in which coupling effects dominate. Although the SFA is not capable of considering the effects of coupling, it allows one to compute the influence of higher order contributions such as quadrupolar modes, which is by definition not possible with the CPDA. Consequently, both approaches provide complementing features that are beneficial depending on the precise use case. Using both approaches also makes it possible to assess whether certain optical features arise due to strong coupling effects or dominant higher order modes in the metasurface.

 

Fig. 27. Tailored Fourier plane with phase randomized image-plane discretization. The simulations were carried out for ${\lambda _0} = 0.6\;{{\unicode{x00B5}{\rm m}}}$ and aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The top row illustrates the specified binary Fourier-plane amplitude distribution with uniformly randomized phase. The metasurfaces are obtained by placing disks at locations where the real parts of the field are positive. The CPDA simulations in the bottom row yield point-symmetric Fourier-plane distributions.

Download Full Size | PDF

H. Tailoring the Angular Scattering

In Section 2.C.1, we have discussed how phase-retrieval algorithms can be applied to solve the phase problem in lensless imaging. Although this is certainly interesting for various applications in the optical regime, this section assesses the applicability of phase-retrieval algorithms for tailored angular radiation patterns. The latter is equivalent to specifying a certain intensity distribution in the Fourier plane and finding the corresponding metasurface that generates this reciprocal intensity distribution. In principle, it is possible to solve this task by Fourier transforming the desired intensity distribution. This is different from lensless imaging, since it is not important to find a unique solution to this problem, and thus the phase distribution in the Fourier plane is not vital. However, it has shown to be beneficial to specify a random phase distribution before Fourier transforming the desired Fourier-plane fields. This is due to the fact that the Fourier transform of an arbitrary Fourier-plane amplitude distribution ${\hat E_0}({q_x}, {q_y})$,

The Fourier transform of the uniformly distributed random phase results in a complex valued random distribution in the image plane. Because of that, the convolution can be understood in terms of a moving average filter of that random distribution with a window function $E({x, y})$, which generally yields a smoother field distribution ${E^\prime}({x, y})$.

A smooth field distribution is often beneficial, if the field distribution has to be converted into a producible metasurface design. Because of that, introducing a random phase in the Fourier plane is recommended as a heuristic principle. A very simple but effective way to convert the obtained field distribution ${E^\prime}(x, y)$ into a metasurface is to threshold the fields such that structures are placed at locations where the real part of the electric field is positive (${\cal R}[{E^\prime}(x, y)] \gt 0$). As illustrated in Fig. 27, only point-symmetric field distributions can be obtained with this simple approach, as the imaginary part of the field distribution is neglected.

A more sophisticated approach to design the angular scattering of metasurfaces is provided by the phase-retrieval algorithms introduced in Section 2.C.1. Since the presented algorithms are iterative, it is usually straightforward to specify further constraints that the metasurfaces must fulfill. For example, the phase-retrieval algorithms intrinsically allow one to target a certain metasurface extent by specifying the size of the support region $S$ accordingly. Zhao et al. [26] have realized a holographic display that can be multiplexed by changing the polarization state of incident light. The presented metasurfaces were designed by modifying the GS algorithm such that the polarization dependence can be incorporated using the Jones calculus. Zheng et al. [70] accomplished sophisticated field distributions with very high diffraction efficiencies and broad bandwidth by using the GS algorithm, but optimizing only for circularly polarized light. The general idea behind employing phase-retrieval algorithms to tailor the Fourier plane is depicted in Fig. 28. Here, it is assumed that an unknown metasurface generates the specified Fourier-plane intensity distribution. Phase-retrieval algorithms can be used to find the image-plane field distribution that corresponds to the desired Fourier-plane pattern. Due to the Shannon–Nyquist sampling theorem, the phase-retrieval algorithm returns a complex valued image with a pixel pitch half the wavelength of incident light. Figure 28(c) reveals that the retrieved image is in very good agreement with the actual metasurface as long as the structural details of the sample are above the diffraction limit. Finally, the metasurface can be generated by distributing custom structures that comply with the amplitude and phase of the retrieved image-plane field distribution at each pixel. Although this has not been shown yet, the CPDA model is a promising tool for improving the design quality of metasurfaces, as it allows one to efficiently compute the Fourier-plane response of extended metasurfaces under consideration of coupling effects, various input polarization states, or even structural perturbations.

 

Fig. 28. Recovery of the image-plane fields via phase-retrieval algorithms. The seed of the pseudorandom number generator to create the positional disorder is zero. The simulations were carried out for ${\lambda _0} = 0.4\;{{\unicode{x00B5}{\rm m}}}$, ${n_x} = {n_y} = 8$, ${p_x} = {p_y} = p$, ${s_{\!\rm{d}}} = 100\%$, ${l_{\rm{c}}} = 400\%$, and aluminum disks with size $0.18 \times 0.18 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. (a) The metasurfaces differ in their initial periodicity $p = {p_x} = {p_y}$ and generate the Fourier-plane amplitudes depicted in (b). The oversampling smoothness (OSS) algorithm retrieves the image-plane field distributions (only amplitudes are shown) (c) from the specified Fourier-plane amplitudes (b) and the support region indicated by the white background in (a).

Download Full Size | PDF

5. MODEL VALIDATION

While the previous section gave an overview on the capabilities of the CPDA model, the goal of this section is to validate the model by comparing the computed results with external software, theoretical predictions, and experiments. Validating the results against three different sources allows one to take advantage of the strengths of each. First, we will compare the spectrally resolved scattering cross section of a single gold disk predicted by the CPDA model with Mie theory. This is particularly useful, as Mie theory provides an analytical solution for a well-defined problem. Subsequently, a slightly more difficult case of coupled gold dimers is considered. Here, the differential scattering cross section ${\sigma _{{\rm{sc}}}}$ at ${\theta _x} = {\theta _y} = 0$ is compared with COMSOL simulations. The primary aim of this step is to assess the normalization and accuracy of the coupling model (see Section 2.A.1) of the CPDA. Then, we validate the image-plane, Fourier-plane, and spectral intensity distributions of largely extended and non-trivial metasurface designs with experimental data. Since no comparable simulational approaches are available that allow one to compute the response of disordered and extended samples with reasonable computational effort, these results cannot be validated using external software.

A. Validation against Mie Theory

Mie theory describes the scattering of electromagnetic waves at spherical objects. The scattering cross section

 

Fig. 29. Validation of the spectral scattering cross sections. (a) The simulations were carried out for a single gold disk with diameter $D$ and shape $D \times D \times 0.03\;{\unicode{x00B5}}{{\rm{m}}^3}$. The normalized scattering cross sections ${C_{{\rm{sc}}}}$ from the CPDA model are in excellent agreement with Mie theory. This indicates that the slight differences in the resonance position of the single particle references in (b) and (c) are due to higher order modes. The simulations in (b) and (c) were carried out for two gold disks with center-to-center distance $d$ and shape $0.05 \times 0.05 \times 0.03\;{\unicode{x00B5}}{{\rm{m}}^3}$. The differential scattering cross section ${\sigma _{{\rm{sc}}}}$ of a single gold disk with the same shape is provided as reference ($d = 0 \;{\rm{nm}}$). All spectra in (b) and (c) are provided for specular reflection at ${\theta _x} = {\theta _y}={ 0^ \circ}$. (b) If incident light is aligned along the dimer axis, the relative amplitudes, resonance widths, and resonance wavelengths agree well with COMSOL simulations. The findings are very similar if the incident light is polarized perpendicular to the structure axis in (c). The dimer differential scattering cross-section spectra are normalized to the amplitude of the single particle spectrum.

Download Full Size | PDF

Strictly speaking, this relation applies only to spherical objects. Nevertheless, we will assume that it also holds true for other shapes such as disks. Thus, we imply that all effects emerging from the altered particle shape are contained within the polarizability.

Figure 29(a) compares the backward scattering cross section of the CPDA model with Mie theory. The scattering cross section of the CPDA model is derived from the Fourier-plane intensity distribution ${I_R}({\theta _x}, {\theta _y})$ and Eq. (76) as

B. Validation against COMSOL

In the next step, we validate the CPDA model against COMSOL simulations for a slightly more difficult system of two gold disks with varying center-to-center distance $d.$ Even so, the differential scattering cross section of a single disk is considered as reference. In contrast to the previous section, the differential scattering cross section ${\sigma _{{\rm{sc}}}}$ instead of the scattering cross section ${C_{{\rm{sc}}}}$ is considered. Both quantities relate to each other as

As indicated by the $y$-labels in Figs. 29(b) and 29(c), only the differential scattering cross section in backwards direction (${\theta _x} = {\theta _y} = \theta ={ 0^ \circ}$) is considered. In Fig. 29(b), the polarization state of incident light is aligned with the dimer. Here, COMSOL simulations reveal that the resonance amplitudes and wavelengths scale highly nonlinearly with the center-to-center distance $d.$ Nevertheless, the CPDA model accomplishes precisely rendering this behavior of line shapes. This is particularly remarkable, as the CPDA model does not consider any extent or shape of the structure, and the edge-to-edge distance for $d = 60 \;{\rm{nm}}$ is only 10 nm. Two minor differences can be found between the models. First of all, the CPDA simulations exhibit sharp features around the resonance positions, which are not present in the COMSOL simulations. It is likely that this is due to numerical errors introduced during the inversion of the transfer matrix ${\hat {\cal T}_M}$ from Eq. (17), because ${\hat {\cal T}_M}$ may become ill conditioned if resonances with small linewidths are present. Second, the single gold disk provided as references is slightly blueshifted in comparison to the COMSOL results. The good agreement of the CPDA approach with Mie theory for the same scenario suggests that the differences are due to higher order modes, which are considered in COMSOL.

 

Fig. 30. Schematic of the experimental setup. (a) The setup used to capture the image-plane appearance consists of an objective with a tube lens and subsequent $4f$ setup that projects the sample onto the RGB CCD camera. (b) To capture the Fourier plane instead of the image plane, the lens ${{\rm{L}}_2}$ is replaced by a Bertrand lens that images the Fourier plane within the $4f$ setup onto the CCD camera. In both setups (a) and (b), the sample is illuminated using Köhler illumination and a white light lamp, which is coupled into the setup with a beam splitter. The $XY$ stage enables to change the incidence angle of the illuminating light on the sample. (c) The metasurfaces are attached to a glass substrate and embedded within the IC1-200 dielectric ($n = 1.41$) to approximate a homogeneous refractive index of $n = 1.5$. An oil-immersion objective with a normalized NA of 0.93, ${n_{{\rm{oil}}}} = 1.518$ and $100 \times$ magnification is used. Spectral data are recorded by replacing the RGB CCD camera with a grating spectrometer. The figure has been adapted from Ref. [16] and depicts the same optical setup. BFP, back focal plane; BS, beam splitter; TL, tube lens; IP, image plane; FP, Fourier plane; BL, Bertrand lens.

Download Full Size | PDF

In Fig. 29(c), the same dimers are shown for perpendicular polarization of the incident beam. The findings are very similar and generally reveal very good qualitative agreement between the CPDA and COMSOL simulations. In this case, less variation of the resonance position is observed. This is due to the excitation polarization, which is perpendicular to the dimer axis and thus weakens the near-field coupling between disks. As before, the CPDA simulations exhibit sharp features around the resonance positions and also the single particle reference simulation is blueshifted in comparison to the COMSOL results.

C. Validation against Experiments

Although any simulation model must ultimately resemble the experimental reality, it is often difficult to use experimental data as ground truth for the validation. This comes from the fact that experiments involve plenty of parts that possess properties that are far from ideal. Many of such properties can be compensated for within the experiment or considered in the simulation model, but there will always remain experimental circumstances that are not included in the model. Such experimental circumstances can be the nonlinear response of the detecting camera, absorption losses and dispersion in the optical components, wavefront errors due to aberrations, and scattering at dust particles in the optical setup, among others. For the experimental setup used in this section, special care had to be taken to embed the metasurface in a constant refractive index environment. Furthermore, matching the color model from Section 2.D to the color model of the CCD RGB camera was a particularly challenging task. Nevertheless, there were no other computational approaches available that allow one to compute the response of plasmonic metasurfaces comprising several thousands of individual structures with reasonable computational effort, so that an experimental validation is the only remaining option.

Figure 30 depicts a sketch of the experimental setup used to acquire the Fourier- and image-plane intensities as well as the spectrally resolved reflectances. Figure 30(c) illustrates that the metasurfaces, attached to a glass substrate, are covered with IC1-200 dielectric and connected to the objective by an immersion oil to mimic a constant refractive index of $n = 1.5$. The objective has a material independent NA of 0.93, so that very high resolutions can be accomplished. Köhler illumination is used for the illumination path to ensure uniform illumination of the sample and enable the independent manipulation of the field and aperture stops. The illumination beam is generated by a spectrally flat laser-driven light source (Energetic EQ99) and coupled into the microscopy setup by a beam splitter located between the objective and tube lens. The imaging path consists of the objective and a tube lens with a subsequent $4f$ system, which projects the intermediate image plane on the RGB CCD camera (Allied Vision Prosilica GC2450) or a grating spectrometer (Princeton Instruments Pixis 256 camera with Princeton Instruments SP-2500i monochromator). The latter depends on whether spatially or spectrally resolved measurements are conducted. The spectral measurements are recorded along an approximately one-dimensional slice in the Fourier or image plane so that spatial information can be captured simultaneously. However, this involves considerably more effort to scan both spatial dimensions. To record the Fourier plane, the ${{\rm{L}}_2}$ lens in the $4f$ setup is replaced by a Bertrand lens, which projects the intermediate Fourier plane onto the RGB CCD camera or grating spectrometer, respectively. For such measurements a zeroth order stop can be assembled in the intermediate Fourier plane, as illustrated in Fig. 30(b).

 

Fig. 31. Experimental validation of the image-plane and Fourier-plane appearances at $0.5\;{{\unicode{x00B5}{\rm m}}}$ periodicity. The seed of the pseudorandom number generator to create the positional disorder is zero. The corresponding metasurfaces consist of ${n_x} = {n_y} = 50$ aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The simulated appearances are derived from the hyperspectral datasets using the color model explained in Section 2.D. The luminosities of the simulations and measurements are matched to render the same visual impression. These results are published in Ref. [16].

Download Full Size | PDF

The fabrication steps of the metasurfaces used for the experimental validation are as follows. First, a positive resist is spin-coated followed by electron beam lithography to define the structure distribution. After developing the resist, the aluminum is deposited onto the substrate with electron beam evaporation. A N-ethyl pyrrolidone (NEP)-based liftoff process removes the residual resist. Eventually, the remaining aluminum nanostructures are spin-coated with IC1-200 to approximate a constant refractive index of $n = 1.5$.

The following description applies to all data depicted in Figs. 31–34. The simulated image- and Fourier-plane appearances are derived from hyperspectral simulations in the wavelength range of $0.38\;{{\unicode{x00B5}{\rm m}}} \le {\lambda _0} \le 0.8\;{{\unicode{x00B5}{\rm m}}}$ so that all perceivable wavelengths are covered. This means that each pixel in the image or Fourier plane has an underlying spectrum, which is converted to the sRGB color space with the color model explained in Section 2.D. The simulation NA is matched with the NA of the objective for better comparability. The white balance of the experimental data is adjusted with a reference aluminum mirror, and all depicted images are background corrected by subtracting an image without a sample from the recorded image of interest. This method of correcting the background is inaccurate but assumed to be sufficient. Both the simulated and experimental appearances have been adjusted in brightness to match the visual impression and bring out all relevant features. The simulated spectra are taken from the zeroth order in the Fourier plane, while the experimental spectra are acquired from reflectance measurements across a central slice of the sample. This difference causes many uncertainties, which is why the spectral data have to be taken with care. An approach that provides better comparability of the spectral data is discussed in Section 6.

 

Fig. 32. Same as in Fig. 31, but for spectral reflectance. It is important to mention that the simulated and experimental data do not show precisely the same spectral quantities.

Download Full Size | PDF

 

Fig. 33. Same as in Fig. 31, but for a perdiodicity of $0.25 \;{{\unicode{x00B5}{\rm m}}}$. These results are partially published in [16].

Download Full Size | PDF

 

Fig. 34. Same as in Fig. 32, but for the spectral reflectance at $0.25\;{{\unicode{x00B5}{\rm m}}}$ periodicity.

Download Full Size | PDF

The first set of measurements was performed for a super-diffraction grating with ${p_x} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$ and ${n_x} = {n_y} = 50$ aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The disorder parameter sets are ${s_{\!\rm{d}}} \in \{0\% , 20\% , 40\% , 80\% \}$ and ${l_{\rm{c}}} \in \{0\% , 200\% , 400\% \}$, whereas only ${l_{\rm{c}}} = 0\%$ is considered for ${s_{\!\rm{d}}} = 0\%$, yielding 10 different realizations.

Figure 31 shows at the top the computationally and experimentally acquired image-plane appearances for different disorder parameters. Generally, the agreement is very good. A slight mismatch in the tone of the simulated appearances can be observed. In comparison to the experimental data, a warmer and more purple tint is present. The differences are particularly present at low disorder strengths. This indicates that the mismatch arises from structural perturbations introduced in the manufacturing process and not considered in the simulation model. These unwanted perturbations get increasingly negligible with growing disorder strength, which leads to better agreement between simulations and measurements. The accuracy of the predicted appearances gets particularly obvious for high disorder strengths and correlation lengths. For these metasurfaces, almost all spatial features are predicted precisely in their extension, hue, and luminosity.

The Fourier-plane appearances are depicted in Fig. 31 at the bottom. Again we find very good agreement between simulations and measurements. The radiation pattern of the simulated and non-disordered metasurface exhibits purple features at the inner edges of the diffraction orders. These features are not present in the measurements and thus suggest that certain spectral regions are overestimated in the simulations. Possible reasons are additional structural defects in the measured metasurfaces, oxidation of the aluminum disks, absorption losses in the optical components, or strong higher order moments not considered in the CPDA model. Besides these purple features, we find that the isotropic scattering that emerges at high disorder strengths is overestimated by the simulations. However, the simulations correctly predict that the (1,0) and (0,1) diffraction orders are still present at ${s_{\!\rm{d}}} = 80\%$ and ${l_{\rm{c}}} = 0\%$, while the (1,1) orders have vanished. Another observation that is true for most of the presented appearances is that weakly pronounced characteristics are less visible in the measurements. A possible reason for this is that the RGB camera can record a maximum of 12 bit, while the simulation model performs the calculations with 64 bit. The higher bit depth allows one to distinguish dark features much longer from an actually black pixel and thus facilitates their visibility.

Before discussing the spectral data shown in Fig. 32, it must be stressed again that the simulations and measurements do not show precisely the same quantities. The simulated intensities are acquired from the intensities in zeroth order of the Fourier plane, while the measurements correspond to the reflectance along a central slice of the metasurface. Thus, perfectly matching results cannot be expected from the outset, and the results must be taken with care. The most obvious difference between the measurements and simulations is that the simulations overestimate the peak that occurs close to ${\lambda _0} = 0.8\;{{\unicode{x00B5}{\rm m}}}$ for low to intermediate disorder strengths. This certainly contributes to the deviations in the predicted appearances to some extent. However, it is questionable how strong the influence of this spectral feature on the appearance actually is, since human perception ends at $0.78\;{{\unicode{x00B5}{\rm m}}}$ and is very weak just below. The measurements generally exhibit a stronger decay towards longer wavelengths in comparison to simulations. Besides that, we often find good agreement between the relative heights and spectral positions of the computed and measured peaks.

The second set of measurements was performed for a sub-diffraction grating with ${p_x} = {p_y} = 0.25\;{{\unicode{x00B5}{\rm m}}}$ and ${n_x} = {n_y} = 50$ aluminum disks with size $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$. The disorder parameter sets are ${s_{\!\rm{d}}} \in \{0\% , 20\% , 40\% , 80\% \}$ and ${l_{\rm{c}}} \in \{0\% , 200\% , 400\% \}$, whereas only ${l_{\rm{c}}} = 0\%$ is considered for ${s_{\!\rm{d}}} = 0\%$, yielding again 10 different realizations.

Figure 33 displays at the top the image-plane appearance of sub-diffraction metasurfaces. The agreement between simulation and measurement is particularly good for the non-disordered case. Both show very similar appearances and exhibit several edge modes with purple and cyan tint. However, effects of a non-uniform illumination, possibly due to pollution of some optical elements, are visible in the measured image-plane appearances, which are additional sources of disturbances not considered in the simulation model. Similar to the super-diffraction metasurfaces, the appearances of the simulated data are shifted towards the purple tones. Furthermore, dark features are again much less pronounced in the measurements, which is presumably caused by the lower bit depth of the RGB camera. Besides that, very good predictive power of the CPDA model is observed. The occurrence of the green areas is accurately reproduced by the simulations. Also, the red and purple domains that appear at high correlation lengths are present in measurements and simulations. However, especially, the purple regions are much less pronounced in the measured data.

Figure 33 displays the Fourier-plane appearance of the sub-diffraction metasurfaces at the bottom. As expected for a sub-diffraction grating, no diffraction modes are visible so that only specular reflection is present in the non-disordered metasurface. For uncorrelated disorder, the simulations predict specular reflections as the dominant feature up to ${s_{\!\rm{d}}} = 40\%$, while at ${s_{\!\rm{d}}} = 80\%$, isotropic scattering sets in. This is different from the measurements, for which specular reflection is dominant for all uncorrelated disorder strengths. However, the simulations accurately predict the desaturation of specular features with increasing uncorrelated disorder. The simulations generally predict a more saturated appearance of the Fourier plane, which comes from differences in the predicted spectral intensities, as can be seen in Fig. 34. Besides that, the simulation model precisely predicts the emergence of a broadened zeroth order for correlated disorder. Both simulations and measurements exhibit slightly more pronounced broadening for smaller correlation lengths and ${s_{\!\rm{d}}} \gt 0\%$.

As explained before, it should not go unmentioned that also the measured and simulated spectral data in Fig. 34 represent slightly different quantities.

Most obviously, the simulations differ from the measurements in the predicted linewidth of the resonance. This is to be expected to some extent, since the manufactured metasurfaces always show imperfections that introduce additional damping losses and thus broaden the resonance. Furthermore, the formation of aluminum oxide within the aluminum nanostructures is commonly known to promote differences between theoretical and experimental data [71]. The narrow linewidth of the simulated data is also the reason for the enhanced saturation in the simulated appearances depicted at the bottom of Fig. 33. Besides that, the resonance position is predicted with good agreement for low disorder strengths. For higher disorder strengths, the measured resonances are redshifted, while the simulated resonances shift to shorter wavelengths. Although the resonance shifts are not predicted correctly, the spectral broadening caused by disorder is reproduced well by the CPDA model.

6. CONCLUSION AND OUTLOOK

This tutorial provided the theoretical basis required to understand and implement the CPDA model. We have presented a disorder formalism that allows one to generate plasmonic metasurfaces with a tunable degree of disorder and correlation. We have shown that the disorder formalism is designed such that the location, orientation, and size of individual nanostructures can be manipulated progressively by setting the disorder parameters ${s_{\!\rm{d}}}$ and ${l_{\rm{c}}}$. Especially, the positional disorder was quantified using the TPCF and SF, which both are established metrics in different fields of research. The SF was connected with the scattering properties of nanostructure ensembles via the so called SFA, which allows one to separate the contributions from the structure distribution and the individual scattering characteristics of the nanostructures. With these disorder metrics, the influence of the different disorder parameters ${s_{\!\rm{d}}}$ and ${l_{\rm{c}}}$ was quantified, and the observed effects were discussed. Since the CPDA is capable of computing spectrally resolved, angular, and spatial far-field distributions, the concepts of the Fourier and image planes were explained using the example of a $4f$ setup. A dedicated color model was introduced that is derived from the sRGB color space and allows one to convert the computed hyperspectral Fourier- and image-plane intensities into color appearances.

A short overview on the implementation of the CPDA package with the associated disorder and color model was given in Section 3. This overview provided insigths on the structure of the developed package and specified the most relevant simulation parameters.

The capabilities of the implemented models were demonstrated, and the consistency of the predictions with theoretical expectations was discussed. It was shown that the model can not only compute the scattering of disordered metasurfaces, but is also able to reproduce the results of recently published metasurface designs [45,72,73] for different polarization states. The CPDA model was compared with the SFA, and excellent agreement was found for the examined metasurfaces. A first step towards truly tailored scattering properties of plasmonic metasurfaces was taken, and the underlying phase-retrieval concepts were explained.

The spectral, spatial, and angular intensities predicted by the CPDA model were validated against Mie theory, FEM simulations, and experiments. While excellent agreement with theory and FEM simulations was found for most of the characteristics, the normalization of the CPDA was shown to be not correct yet. The agreement between the CPDA and experimental data of the Fourier- and image-plane appearances was mostly very good and exhibited only a slight mismatch in the global tint. Dark features, such as isotropic scattering backgrounds, were generally less pronounced in the experimental data. The spectral agreement with the experiment was satisfactory but showed deviations in the predicted linewidths and amplitudes. However, many of these deviations were traced back to differences in the data acquisition, fabrication inaccuracies, and oxidation of the aluminum nanostructures.

The CPDA model has demonstrated great predictive power for a variety of disordered and very sophisticated plasmonic metasurfaces. As shown in Fig. 35, it is capable of simulating the monochrome scattering properties of metasurfaces with 10,000 nanostructures in only a few minutes with good resolution. The run time is mainly determined by a quadratic dependence on the number of structures [see Eq. (17)] and a linear dependence on the number of pixels and wavelengths. The memory requirements are also governed by the quadratic dependence on the number of structures to store the transfer matrix ${\hat {\cal T}_M}$ and are mostly independent of the size of the simulation domain. In this tutorial, most simulations were run on a Windows machine with 16 GB RAM and four cores. The simulations in Figs. 31–34 were computed on a Linux machine with 128 GB RAM and 16 cores for speed reasons. However, assuming that the transfer matrix consists of $m = 128 \;{\rm{bit}}$ double precision complex entries for ${n_{\rm{s}}}$ structures in ${n_{{\dim}}} = 3$ dimensions, it is possible to roughly estimate that on a nowadays common $M = 16 \;{\rm{GB}}$ desktop machine, simulations with up to ${n_s} = 10{,}000$ structures should be possible [$M = m \cdot {({n_{{\dim}}} \cdot {n_{\rm{s}}})^2}$]. This makes it a very powerful tool for tackling the challenge of tailored scattering under the consideration of structural defects and coupling effects present in real systems. Due to the integration of the sRGB color model, the CPDA is not only applicable for physical computations, but additionaly enables to derive the spatial and angular appearances of plasmonic metasurfaces.

 

Fig. 35. Computation times for varying numbers of structures and grid resolutions. The mean and standard deviations of the image- and Fourier-plane computations at ${\lambda _0} = 0.6\;{{\unicode{x00B5}{\rm m}}}$ are displayed. The resolution of the image- and Fourier-plane grid is ${n_{{\rm{px}}}} = {n_{{\rm{px}},x}} \cdot {n_{{\rm{px}},y}}$. Each set of simulation parameters is computed five times to derive the statistical parameters. All simulations were carried out for ${n_{\rm{s}}}$ aluminum disks of shape $0.12 \times 0.12 \times 0.06\;{\unicode{x00B5}}{{\rm{m}}^3}$ and randomly distributed over intervals ${-}{p_x} \cdot ({n_x} - 1)/2 \le x \lt {p_x} \cdot ({n_x} - 1)/2$ and ${-}{p_y} \cdot ({n_y} - 1)/2 \le y \lt {p_y} \cdot ({n_y} - 1)/2$ for ${p_y} = {p_y} = 0.5\;{{\unicode{x00B5}{\rm m}}}$ and ${n_x} = {n_y} = \sqrt {{n_s}}$. All quadratic polynomial fits (solid lines) have ${R^2} \gt 0.999$ and thus confirm the quadratic dependency on the number of structures. These simulations are executed on a Linux machine with 128 GB RAM and 16 cores to reach high values of ${n_{\rm{s}}}$. However, note that most simulations here have been calculated on a standard desktop.

Download Full Size | PDF

Nevertheless, many tasks remain to be done in future work. Some of these tasks revolve around the completion of the CPDA, while others involve the application of the approach on new problems. First and foremost, the correct normalization of the model has to be found so that all simulated properties can be related precisely to the physical reality. Also, the reciprocity theorem can be formulated such that it provides a solution for the forward scattering (transmission) of plasmonic metasurfaces. This is interesting not only as it allows one to predict the transmission characteristics, but also because the absorption within the metasurface can be deduced from the computed forward and backward scattering. Another measure that would allow one to compute a wider range of metasurfaces is to extend the database of available nanostructure polarizabilities. This particularly concerns the polarizabilities of bigger structures, since a large body of published metasurfaces relies on those. Besides these rather straightforward steps, there are also some steps that are certainly more difficult to implement but still very interesting. First of all, the model would tremendously benefit from the implementation of a substrate beneath the metasurface for the simple reason that the substrate is indispensable in experimental systems. Generally, a formalism that allows one to efficiently incorporate layered structures would yield a model with much broader applicability.

Another yet untackled task is to consider magnetic modes in the CPDA model. This is required to accurately predict the scattering of dielectric metasurfaces, which often exhibit strong magnetic modes. Besides these issues directly related to the CPDA, the quantification of the rotational and dimensional disorder with suitable metrics remains to be done.

Although the results have not been validated at this point, it is possible to examine three-dimensional structure distributions with the implemented CPDA. However, even if three-dimensional metamaterials would certainly yield interesting scattering distributions, it is questionable whether those designs can be fabricated with reasonable effort. From an experimental perspective, the application of the CPDA model on structure distributions with unique long- and short-range orders is of great interest. Hyperuniformity [67–69], Matérn type processes [74,75], or random sequential adsorption [76] are very promising disorder-based concepts for enhancing the conversion efficiency of solar cells and creating flexible interfaces, among others [19]. While these structure distributions come from theoretical considerations about order, very sophisticated scattering distributions can be obtained from the application of phase-retrieval algorithms [26]. The underlying concepts have been discussed in this tutorial, but the proof that the CPDA model is a very effective tool for the precise design of angular scattering is still pending.

There is also space for improvement regarding the experimental validation of the CPDA model. Most importantly, the validation of the appearance should not be carried out anymore by directly measuring the RGB appearance. Instead, it is much more reasonable to measure the monochromatic Fourier- and image-plane intensity distributions and derive the appearances with the same color model used for the simulated data. This way, it is not necessary to consider any differences caused by deviations in the illuminating light source, inaccurate white balances, or differences in the employed color models. Additionally, it allows one to directly compare the measured and simulated scattering characteristics at different wavelengths. Furthermore, the experimentally acquired spectral data should be measured in the zeroth order of the Fourier plane. This ensures that there are no avoidable ambiguities in the interpretation of the experimental data. It is also advisable to use gold structures in prospective experimental validation measurements, since this avoids the problem of oxidation that aluminum suffers from.

So far, the fabricated metasurfaces are derived from metasurface designs passed to the CPDA model. However, this comes with the disadvantage that inaccuracies originating from the manufacturing process cannot be considered in the simulations. It is conceivable to overcome this problem by extracting the precise location, orientation, and size of the fabricated metasurfaces from scanning electron microscope (SEM) images and feeding these into the simulation model.