3D Plasmonic Design Approach for Efficient Transmissive Huygens Metasurfaces

In this paper we present a design concept for 3D plasmonic scatterers as highefficiency transmissive metasurface (MS) building blocks. A genetic algorithm (GA) routine partitions the faces of the walls inside an open cavity into aM x N grid of voxels which can be either covered with metal or left bare, and optimizes the distribution of metal coverage needed to generate electric and magnetic modes of equal strength with a targeted phase delay (Φt) at the design wavelength. Even though the electric and magnetic modes can be more complicated than typical low order modes, with their spectral overlap and equal strengths, they act as a Huygens source, with the accompanying low reflection magnitude. Square/hexagonal voxels inside square/rectangular cavities are thoroughly analyzed for operation at 8 μm, although the technique can be applied to different cavity geometries for operation across the electromagnetic spectrum. Results from full-wave simulations show the GA routine can repeatedly pinpoint scatterer geometries emitting at anyΦt value across 2π phase space with transmittances of at least 60%, making these MS building blocks an attractive plasmonic alternative for practical optical applications. Full-scale metasurface devices are calculated from near-fields of the individual elements to validate the optical functionality. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
The engineering of a surface for localized control of the phase of emitted light has been employed for decades; known early developments resided in the radio frequency (RF) regime for the advancement of low-mass beam control [1]. In the past seven years, however, these concepts have re-emerged as a topic of interest under the overarching "metasurfaces" moniker, with simulation and fabrication advances encouraging developments at frequencies previously considered infeasible. Since the seminal two-dimensional (2D) plasmonic works in infrared (IR) and telecom wavelengths [2][3][4], metasurface research has expanded to nearly every possible corner of electromagnetic design. In terms of spectrum, efforts have spanned RF [5][6][7], THz [8,9], IR [2,10,11] and near-IR/visible regimes [12][13][14]; in terms of application, aside from the multitude of lensing and beamsteering functions, metasurfaces have been implemented for other novel wavefront manipulations, such as compound lensing, optical vortex beam formation, and polarization conversion [15,16]; in terms of constituent material base, plasmonics have generally ceded popularity to dielectrics, as interests have shifted towards optical applications in the visible regime [17]; and finally, in terms of architecture, 2D works dominated early, whereas high-aspect ratio dielectrics [18][19][20] and stacked multi-layer pseudo-3D plasmonics now reign.
Despite the wide range of materials, architectures and designs, a primary challenge still remains in attaining an efficient comparison to conventional bulk optics. Single layer 2D plasmonics have fared the worst in this, as it can be shown by applying the surface equivalence principle to Fresnel coefficient equations that mono-modal scatterers illuminated by linearly polarized light cannot attain the full 2π phase control necessary for beam control, nor can any multi-modal scatterer reach an efficiency greater than 25% if linear-to-linear polarization conversion is required [21,22]. On the other hand, much higher efficiencies have been reported for multi-layer plasmonic structures, with the most successful of these having incorporated multiple Huygens-like source elements. A Huygens source is derived from the equivalence principle, which states that a surface exhibiting equal electric and magnetic current amplitudes with appropriate phase can completely eliminate backward wave propagation [23]. These behaviors are emulated in a metasurface through the subwavelength superposition of electric and magnetic responses, demonstrated in various designs using independent capacitive/inductive circuit elements or cascaded sheet admittances acting as control mechanisms [8,[23][24][25][26][27][28][29][30][31][32][33][34].
However, a majority of Huygens metasurfaces have been designed for operation at THz and lower frequencies. At higher frequencies, Ohmic loss is significant and often the lumped circuit elements utilized are dimensionally challenging to fabricate or the precise control of resistive/reactive elements are functionally difficult to implement (e.g., placing passive resistive elements into LC circuits to control phase [27]). Due to these challenges, attempts to attain full phase coverage at IR frequencies using plasmonic Huygens' sources have been sparse, and have often required a different design approach: in addition to attempts at the sheet admittance approach [32,33], optical lumped LC circuits [20], nanoparticle arrays [33] and coupled waveguides [35] have been proposed to provide the out-of-plane (OOP) features necessary for proper mode generation. None of these metallic Huygens source designs have delivered a high transmittance (> 40%) across the entire 2π phase space, in practice; therefore, a continued pursuit of metasurface design concepts is warranted.
In this work, we present an alternative design approach for a plasmonic metasurface architecture, one which does not need explicit modal controls input by a designer and is viable at high frequencies. The metasurface design is focused on a "blank slate" approach, where there are no a priori constraints on the manner in which the electric and magnetic modes are formed or how many modes are necessary to produce the desired phase retardation. This is executed via a 3D unit cell that is populated with voxelated binary grids both in-plane and out-of-plane, and a genetic algorithm (GA) routine which evolves the voxel material into the scatterer topologies needed for generation of electric and magnetic modes to provide the desired performance metrics (target phase retardation, maximum transmittance, and minimum reflectance). This means multiple modes of each kind can coexist locally within a single unit cell, and the result is a metasurface element that looks nothing like a conventional scatterer, topologically. The concept that a GA implementation can develop effective-but exotic-scattering elements has been proven, notably in works on 2D planar metamaterials [36][37][38][39][40][41]; however, this is the first known implementation of this theory for a realizable 3D metasurface element, and boasts an extremely modular design, extendable to an arbitrary number of grids made of arbitrary voxel tessellations inside an arbitrary volumetric geometry.
The transmissive elements can attain any phase value by design, spanning the entirety of the 0-2π phase space and with simulated transmission efficiencies of at least 60%, and thus have great potential to serve as building blocks for metasurface applications. While our examples are presented at long-wave infrared (LWIR) wavelengths, the approach is scalable to any spectral regime. We first examine the proposed design methodology; then, we validate the GA routine numerically through an example design optimized for a 180• phase shift from the 0• phase reference and demonstrate the capability of the GA routine to target any arbitrary phase.

Unit cell design
While the concept of a plasmonic Huygens source is straightforward (e.g., see supporting information for [27,33]), the means of modifying a specific architecture to access a full 2π phase space can be challenging. The solution presented in previously mentioned works is generally the same conventional framework: control the capacitive/electrical response via dipole-like surface currents along the sheets, and control inductance/magnetic response via circulating currents between (planar) or along (out-of-plane) the sheets. Design methods which do not begin with a pre-supposed solution can also be used to design the meta-atoms. Neural networks can be used to design atoms [42]. Rather than assuming some determinable combination of dipole lengths and loop radii are ideal for phasing the Huygens' sources, we decided to pursue an evolutionary optimization approach, where the scatterer begins as a random geometry and evolves into the most effective design by means of a genetic algorithm, in a binary manner similar to previous electromagnetics approaches [36][37][38][39][40][41][42].
The unit cell model employed for this endeavor is based on a 3D fabrication phenomenology termed "membrane projection lithography" (MPL), developed at Sandia National Laboratories [43,44]. The baseline structure consists of an array of silicon (Si) boxes with cubic cavities of air open on one side. In principle, up to four vertical internal faces and the horizontal floor can be decorated with gold (Au) scatterers. Each face is broken into a grid of voxels with a given resolution, with a single voxel being made of either air or Au and representing the fundamental building block for constructing the scatterer geometries. To show the available diversity in the design architecture, Figs. 1(a), 1(b) and 1(c) give examples of a 7 × 9 square grid with rectangular spacing, a 7 × 10 hexagonal ("hex") grid containing n = 82 voxels (counting half-voxels) inside a cubic cavity and a 7 × 9 offset square ("brick") grid containing n = 68 voxels inside a rectangular cavity, respectively. Both the hex grid and brick grid avoid the issue of having "corner-to-corner" voxels touching. Fabrication resolution limits prohibit the formation of sharp corners, making it difficult to resolve such features, so that designs that depend on these features can suffer by the formation of bridges or voids in the as-fabricated structure. By opting for the hex and brick patterns, these patterns are avoided.
Each face in the interior of the unit cell could be independently optimized, potentially resulting in up to five iterations of the MPL fabrication process, a time-consuming endeavor. One approach to manage the fabrication complexity is to use designs which require a single membrane pattern, evaporated multiple times within the unit cell. For instance, in Fig. 1(a) the two vertical faces and the floor are decorated with a single pattern as a demonstration. Note that the pattern structure is projected, so that for the right and left faces the top and bottom of the pattern are reversed. In this case, even though three faces are decorated, only a single pattern is optimized. The simplest alternative is to optimize a single pattern for a single face. For the balance of this work, only a single face will be decorated with plasmonic structures.
Since it was expected that this MPL-based 3D design approach would generate non-analytic solutions with complex geometries, the finite element analysis (FEA) software suite COMSOL Multiphysics with its LiveLink for MATLAB module was employed to generate the total field phase and amplitudes needed for the MATLAB-based GA calculations. Thus, a specific design is evolved through the interplay of the FEA and GA processes, as summarized in the following. A large initial population of models ("individuals") are generated, each with random allocations of voxels ("chromosomes") on any or all faces. Each individual is uniquely defined by a string of n parameters, with each element representing a single voxel and possessing a binary value of either "1" for Au or "0" for air. The initial population is evaluated, performing an infinitely periodic unit cell analysis to extract the S-parameters for the total field. Their transmittance (T = |S 21 | 2 , co-polarized transmission) and phase values (Φ = ∠S 21 ) are compared against a two-valued fitness function which seeks to maximize T at some targeted phase value (Φ t ), given as: where w is the weight applied to Φ or T, σ is a standard deviation defining a Lorentzian distribution centered at Φ t , Φ 0 /T 0 is the phase/transmittance of the bare (unmetalized) Si cavity array, and T low is a minimum acceptable value for T, chosen to be 0.3. Additional constraints were applied to Eq.(1), setting F = 0 for any solutions to mitigate the possibility of solutions with high T /low Φ and low T /high Φ from dominating. This accounts for the null region T > 0.7 As can be seen from an example modified fitness function in Fig. 1(d) using the values Φ t = 120• (red, dashed), Φ 0 = 62• (black, dashed), solutions of constant F form contours (black, solid) of increased curvature towards the desired solution (F → 1). The solution trajectory across generations follows regions of constant fitness (denoted by the '=') and regions where it cuts across contours toward increasing fitness (+). The population is ranked according to highest fitness F and competes in a steady-state, tournament-style selection process, where ideally the fittest individuals are chosen as parents for creation of a new generation of children. A single-point crossover takes a fraction f of a string from one parent and a fraction (1f ) from the other parent to form the new child string; this process is repeated N pop times to form the new generation, and then a mutation is applied to the child string, randomly flipping a low percentage of the bits (nominally 2-3%) from 1→0 or 0→1 . The result is a new set of voxel allocations which are fed into the FEA computation as new models and solved, iterated N iter times until a threshold variance from Φ t and/or T 0 is reached.

Validation of genetic algorithm and MPL design
To validate the binary GA optimization, a series of elements were produced that enable full phase coverage at a high transmittance. Using the reference phase of the unmetalized MPL cavity of Φ 0 = 62•, the 2D solution spaces of a 8 hex layout in a rectangular cavity for an 8-element discretization of Φ t are arrayed in Figs. 2(a)-2(f), showing relative phase shifts from Φ 0 of ∆Φ =(a)-150°and -105°, (b) -60°, (c) -15°, (d) 30°, 75°and 120°, and (e) 165°. Element numbers are linked to their grid layout below (blue = metal), with elements 1, 2 and elements 5, 6, 7 taken from the same GA runs, as shown in Figs. 2(a) and 2(d). Dashed lines show the baseline transmittance of the unmetallized cavity T 0 = 0.746 (red), Φ 0 (blue) and Φ t (black), while yellow dots indicate the performance metrics of the optimal elements chosen for experimental validation, presented in a linear plot in (f). The cross-polarized transmission (T cp )is negligible, indicated by the dashed red curve near the origin in Fig. 2. (f). The grey dots in Fig. 2. (a)-Fig. 2. (e) are intermediary generations. All of the GA runs have very similar characteristics: the same clustering of initial populations around Φ 0 and T 0 and final populations (white dots) around Φ t ; large groupings from which a near-continuum of Φ values can be extracted (improving computation efficiency), and formation of a loose Pareto front-keep in mind some fronts are less well-defined than others due to hitting a termination threshold before it could solidify. Generations will evolve towards + 180°or -180°, depending on which direction is the smallest ∆Φ; as can be seen in (e), this sometimes means passing over the negative phase boundary to reach a positive phase value, and vice-versa. The lowest transmittance was < T> = 0.60 (#7), while the average transmittance was < T> = 0.66 and the average reflectance was < R> = 0.03. Finally, in this work we did not enforce fabrication-specific constraints. Some of the evolved designs had closed loops in them which cannot be suspended during fabrication. These designs are the dots on the target phase line with transmittance greater than the chosen design, and were eliminated by inspection. Ultimately, these fabrication constraints can be built into the evolution process. The source of the phase shifts can be explained through examination of the near-fields around an example 7 × 10 hex scatterer targeted at Φ t = 65°, shown in Fig. 3 the (a) H x and the (c) E z field components around the metasurface interface. Insets provide 4X scaling of (b) H x and (d) E for an enlarged unit cell, along with current densities (arrows). Wherever there is strong rotation in the induced surface currents, there exists strong magnetic modes which can exhibit strong spatial variation in phase; likewise at narrow gaps in the scatterer geometry there are strong electric modes. By the relation Z = E z /H x for a propagating plane wave, we can scale the H-field by the impedance Z in terms of V /m, and we find that the mode strengths of the two components are of the same order. With the close proximity of the electric and magnetic modes, and their near-equal strengths, the Huygens criteria are achieved, and the resulting transmitted total field (left of the unit cell in (a),(c)) is also of near-equal strength. These factors, along with the demonstrated high transmittance and low reflectance (T = 0.63/R = 0.02 for this example), strongly indicate the elements are behaving as Huygens' sources. As seen in Figs. 2(a)-2(e), early generations tend to swarm around the baseline solutions, as sparsely populated grids will tend to poorly exhibit the resonant coupling needed for phase retardation. However, as the voxels begin to evolve into capacitive and inductive geometries, coupling improves and the continuum of Φ/T values becomes more accessible. This can be seen by tracking the evolution of the individual with the highest F over several generations, as seen in Fig. 4. The design layout (gold) of the randomized initial (i = 1) population is the least ordered, and provides very few geometries for interaction with the incident light. As the grid evolves to more contiguous topologies, with out-of-plane curvatures and close-packed capacitive sub-elements, the localized modes appear. The GA routine advances designs where these modes meet the phase and amplitude requirements, and suppress designs which do not. The parameter space is vast: in this case, the hex layout consists of 82 unique voxels, and each voxel is a parameter with a possible value of 1 or 0. To sweep the entire space of possible arrangements would, by the Rule of Product, require 2 n ≈ 5 × 1024 solutions! The GA routine generally solves in 2000-4000 individuals, or about 30-40 generations. While not an exhaustive search of the space, given these considerations, our approach is an efficient means to generate acceptable non-intuitive solutions out of an enormous solution space for a complex, 3D design in relatively few iterations.

Full-scale device simulation
Full-scale metasurface devices cannot be efficiently modeled from such complex individual unit cells using FEA solvers, so to simulate the response of a large array a vectoral diffraction Fig. 4. Evolutions of the voxel layouts (gold) from the highest F designs from an example GA routine, from the initial population (i = 1) to where the "best" solution was determined at generation i = 39. The transverse H x field strength is underlaid. As voxel density increases and coalesces into the early framework of scattering geometries, mode formation and strength increases, and greater phase shifts are attained.
formulation of the Huygens principle was employed, namely a Stratton-Chu far-field projection: where the tangential fields' are taken from a homogeneous, closed boundary S, a distance λ 0 /4 away from the metasurface,n is the unit normal to S, r , is the vector from the origin to S, Z is the impedance µ/E, andr is the unit vector from the origin to the point r. To ease the calculation, the surface integrals can be discretized into i j uniformly-spaced (∆x, ∆z,) summations. For projection onto a 2D-plane at x = 0 from a 1D metasurface, the field contribution from the k th element in the device located at z k can be calculated at any point in the plane: The resulting intensity profiles of (a) a 2 mm beamsteerer and (b) a 3.8 mm diameter lens are shown in Fig. 5, operating at 8 µm. The beamsteerer consists of 108 supercells of eight unit cells based on a square voxel tessellation and designed to emit at 23.5°to the surface normal, with each supercell linearly arrayed in pairs of four distinct phase discretizations spanning 0°, 90°, 180°, and 270°relative phase shifts, whereas the lens was based on an 8-element hex voxel tessellation and spans 3 complete phase cycles of 2π, with a designed focal length of f 0 = 7.5 cm. While the beamsteerer emits at 25.8°(9% error) and the lens focuses at 7.35 µm (2% error), both discrepancies can be attributed to the truncation of the device, and show trends of improvement as the diameter increases. The lens is diffraction-limited to within 0.1%.

Conclusion
In this paper, we proposed novel 3D plasmonic scatterer design concepts for use as building blocks in practical metasurface applications. Using the real-world MPL structure as a basis, we showed that by mating a GA optimization routine into a full-wave FEA solver we could generate a scatterer with a Huygens-like response across the entire 2π phase space, and calculated full-scale devices to confirm proper optical function. When normalized to the scattering losses of the unmetalized rectangular cavities, transmittances of at least 60% were attained for any targeted phase value. This is an unprecedented result for transmissive plasmonic metasurfaces operating at high frequencies, and can be attributed to the OOP voxel distributions allowing multiple distinct, but spatially-localized modes operating at near-equal strengths, disrupting the notion that singular electric and magnetic modes must be controlled through independent mechanisms. Our proposed OOP metasurface building blocks have great potential to produce highly-efficient and extremely versatile low-profile optical devices.