## Abstract

Spatial transforms are a popular technique for designing periodic structures that are macroscopically inhomogeneous. The structures are often required to be anisotropic, provide a magnetic response, and to have extreme values for the constitutive parameters in Maxwell's equations. Metamaterials and photonic crystals are capable of providing these, although sometimes only approximately. The problem still remains about how to generate the geometry of the final lattice when it is functionally graded, or spatially varied. This paper describes a simple numerical technique to spatially vary any periodic structure while minimizing deformations to the unit cells that would weaken or destroy the electromagnetic properties. New developments in this algorithm are disclosed that increase efficiency, improve the quality of the lattices and provide the ability to design aplanatic metasurfaces. The ability to spatially vary a lattice in this manner enables new design paradigms that are not possible using spatial transforms, three of which are discussed here. First, spatially variant self-collimating photonic crystals are shown to flow unguided waves around very tight bends using ordinary materials with low refractive index. Second, multi-mode waveguides in spatially variant band gap materials are shown to guide waves around bends without mixing power between the modes. Third, spatially variant anisotropic materials are shown to sculpt the near-field around electric components. This can be used to improve electromagnetic compatibility between components in close proximity.

## 1. Introduction

Electromagnetic fields cannot be controlled inside a single homogeneous material. Real devices must contain gradients, curves, interfaces or other inhomogeneities that interact with the field to produce a response. In this regard, photonic crystals [1,2] and metamaterials [3,4] that are macroscopically homogeneous are very limited in what they can be made to do. Transformation optics (TO) [5–7] has quickly become the dominant technique for designing periodic structures that are macroscopically inhomogeneous. The technique takes a spatial transform as the input, applies it to Maxwell's equations, pulls the transform out of the spatial coordinates and incorporates it into the constitutive parameters. The output of the technique is a map of permeability and permittivity as a function of position. The permeability and permittivity functions are typically tensor quantities containing extreme values that only metamaterials and photonic crystals are able to provide. In order to adjust the effective constitutive values throughout a lattice, it is necessary to spatially vary, or functionally grade, different attributes of the lattice, such as the structure within the unit cell, fill fraction, lattice spacing, orientation of the unit cells, material composition and even the symmetry of the lattice. Generating the geometry of spatially variant lattices, however, is challenging because the electromagnetic properties of the lattice depend heavily on the size and shape of the unit cells and often require adjacent unit cells to be the same. This means that lattices must be twisted, bent and reshaped while keeping each unit cell consistent with the design. Further, this must be accomplished without distorting the unit cells and without introducing unintended defects or their electromagnetic properties will be weakened or destroyed altogether. The ability to keep the geometry of the unit cells in a spatially variant lattice consistent with that of a uniform lattice is key to maintaining the same electromagnetic properties as the uniform lattice. To date, this type of spatial variance has only been accomplished in simple and canonical configurations. The vast majority of this work spatially varied only the fill fraction [8–21] because this is the simplest to do from a mathematical perspective. Other work can be found that spatially varied the lattice spacing [22,23], the orientation of the unit cells [24–27], the pattern within the unit cell [28–30], the symmetry of the lattice [31] and combinations of these [32–35]. The concept of spatially variant lattices applies equally well to metasurfaces and work can be found in the areas of leaky-wave antennas [36], radiation and guiding [37], phase gradient metasurfaces [38,39], Luneberg lenses [40–42] and more [38,39,43].

A four-step design process using spatial transforms is summarized in figure 1. The first step is to define a coordinate transform that describes how the fields or waves should be bent or directed in space. The second step applies the spatial transform to Maxwell's equations, but moves the transform out of the spatial coordinates and into the constitutive parameters *μ* and *ε*. This step is what is called TO and it produces a map of the permeability and permittivity as a function of position. It says nothing about how those constitutive parameters are to be realized. This information comes from step 3, where the constitutive values generated by TO are linked to engineered materials that can provide the required material properties, even if only approximately. Each point is linked to a potentially different metamaterial, different orientation, different lattice spacing and even a different lattice symmetry. In this sense, step 3 converts the maps of *μ* and *ε* into maps of geometrical parameters that are to be spatially varied throughout the device. Given this information, the final step is to generate the overall spatially variant lattice that is smooth, continuous, free of defects and that correctly realizes the required geometry at every location in the lattice.

The focus of this paper is on the algorithm associated with step 4 in figure 1 and on new device concepts that the algorithm enables. To date, the algorithm described by Rumpf & Pazos [44] is the only general-purpose procedure presented in the literature for generating the physical geometry of spatially variant lattices. Appendix A summarizes this method and presents some new developments for generating better quality lattices and for generating them more efficiently. The method is presented here in the context of three-dimensional lattices, but the algorithm is equally well applied to two-dimensional metasurfaces. For example, in prior work, a guided-mode resonance filter was compensated to work on a curved surface by spatially varying the period of the grating [45]. The ability to spatially vary unit cell orientation, spacing and pattern throughout a periodic structure enables many new design possibilities. In addition to being a necessary design step when using spatial transforms, the ability to generate high-quality spatially variant lattices enables entirely new design paradigms that cannot be conceived using spatial transforms. When the TO step is not used, the input to step 4 in figure 1 must be constructed by other means. Often this entails constructing a map of the Poynting vector to direct the flow of waves along some path. Three new design paradigms are discussed in the following sections. First, spatially variant self-collimating photonic crystals were demonstrated to flow unguided beams around a 90° bend more abruptly than anything else reported in the literature while also using a much lower refractive index. New results are presented to realize smaller lattices operating at shorter wavelengths. Second, simulation results are presented for multi-mode waveguides in spatially variant band gap materials that guide waves around sharp bends without mixing power between the modes, thus maintaining the modal configuration in the waveguide. Third, spatially variant anisotropic metamaterials (SVAMs) are shown to be able to sculpt the near-field around devices. This can be used to improve electromagnetic compatibility (EMC) of multiple components placed in close proximity.

## 2. Spatially variant photonic crystals

At a macroscopic level, devices designed by TO are based solely on refraction, so they are limited in this regard. Periodic structures provide mechanisms besides refraction that can provide more abrupt control of electromagnetic waves using more realistic material properties. This section discusses two concepts using spatially variant photonic crystals (SVPCs) that outperform any device to date designed using spatial transforms. Photonic crystals, however, are based on resonance so they are fundamentally bandwidth limited. In principle, devices designed with spatial transforms are not bandwidth limited because they are based on refraction, which is not sensitive to frequency in the absence of dispersion. However, the underlying periodic structures producing the effective properties often rely on resonant phenomena that are bandwidth limited.

### (a) Abrupt bends using spatially variant self-collimation

Self-collimation is a phenomenon observed in photonic crystals where beams propagate without diverging and are forced to propagate along an axis of the lattice [46–51]. The phenomenon has also been applied to surface waves [52] including spoof surface plasmons [53,54]. It follows from the theory that beams can be made to flow along arbitrary paths inside these lattices by spatially varying the orientation of the unit cells [25,55,56]. These lattices cannot be designed by TO because the effective material properties are constant throughout. It is only the orientation of the unit cells that is adjusted. Self-collimation is based on diffraction instead of refraction so it can control waves over more abrupt length scales without requiring extreme material properties. To construct the input maps for the spatially variant lattice synthesis algorithm, it was only necessary to define the Poynting vector as a function of position. The map for a 90° bend is shown in the leftmost diagram of step 1 in figure 1. In this case, it was only necessary to spatially vary the orientation of the unit cells, not the lattice spacing or fill fraction. No coordinate transform was needed and TO was circumvented in this manner.

The first demonstration of a self-collimating SVPC was at microwave frequencies [56]. This device, shown in figure 2*a*, operated at around 15 GHz and flowed an unguided beam around a sharp 90° bend. Soon after, an SVPC with lattice spacing of 2.0 μm was demonstrated at optical frequencies at a vacuum wavelength of 2.94 μm [55]. This device is shown in figure 2*b*. The bend radius for both devices was around 6.4λ_{0}. The three-dimensional nanophotonic lattice was fabricated by multi-photon direct laser writing (DLW) in the photopolymer SU-8. DLW is a nanoscale digital manufacturing technique similar to three-dimensional printing, which can be used to create free-standing three-dimensional structures in a variety of materials, including common photopolymers. The method of DLW has been thoroughly described in [57–59] and the details of our particular implementation are described in [60]. At near-infrared wavelengths, SU-8 has a refractive index of around 1.57 [55]. This is similar to that of other photopolymers, but it is small compared with the refractive index most typically used for photonic crystals. Through simulation, it was determined that the best performance was expected for a volumetric fill fraction of roughly 50% [49]. This was later confirmed experimentally.

Figure 2*b* shows a scanning electron microscope image of an optical SVPC with a volumetric fill fraction of 52%. The SVPC is viewed from the top looking down onto the lattice and onto the supporting substrate. Its optical performance was measured by scanning optical fibres along the three exit faces of the SVPC, while a fourth fibre introduced source light onto the input face. The fibres were coupled to infrared detectors so the intensity of the light could be measured. The blue lines in figure 2 give the measured irradiance relative to the input beam as the receive fibres were scanned across the SVPC faces. These data show that the SVPC effectively flows light towards the ‘bent’ face. The peak-to-peak ratio of the light exiting the ‘bent’ face to the light exiting the ‘straight’ face was 8.5 : 1. For this particular SVPC, a strong signal is also seen at the ‘opposite’ face. From simulations, it was determined that this was caused by scattering at the edges of the lattice. More work is needed to suppress this scattering and a number of techniques have already been described in the literature [11,61,62]. Overall, very little light is passing straight through the SVPC, which indicates that this device can flow light through tight turns as designed.

The beam-bending efficiency can be optimized by adjusting the volumetric fill fraction of the SVPC. Figure 3 shows three different SVPCs fabricated with different fill fractions along with the resulting irradiance line scans measured at the ‘bent’ face. The beam-bending efficiency can be gauged from the peak-to-peak ratio of the irradiance scans measured at the ‘bent’ and ‘straight’-through faces. Below a fill fraction of 40%, the peak-to-peak ratio is low and just above unity. At fill fractions above 55%, the peak-to-peak ratio is also low, falling below 0.5%. The highest beam-bending efficiency, with a ratio exceeding 8.5 : 1, is achieved for the SVPC having a fill fraction near 50%. Other functional aspects such as the angle through which the beam is bent, the degree of power splitting and the polarization selectivity can all be controlled through the design of the SVPC [55].

Figure 4 shows two SVPCs that have an axial feature width of only 320 nm and axial unit cell spacing of 720 nm. These were fabricated by DLW in the photopolymer IP-Dip (Nanoscribe) using careful control of the exposure conditions to achieve the small size. This structure will bend light having λ_{0}=1.0 μm. Further improvements in the DLW processing and in the performance of the photopolymers should enable fabrication of SVPCs that can control light at λ_{0}=850 nm. This is the output wavelength of an important class of vertical cavity surface-emitting lasers (VCSELs). Controlling the output of these VCSELs is of interest for enabling new architectures for chip-to-chip optical interconnects and optical networking [63].

Controlling deformations becomes critical when a beam is to be flowed around more than one bend because the spatial variance is more complicated. To demonstrate this, two SVPCs were generated using the self-collimating lattice described in [64] to flow a beam through two very tight consecutive 90° bends. These were simulated using the finite-difference frequency-domain method [65], which has proved to be a highly effective tool for simulating devices incorporating spatially variant lattices. The SVPC shown in figure 5*a* was generated using the basic algorithm. The self-collimation effect was weakened by deformations that were most severe near the centre of the lattice. This caused scattering and reduced transmission to below 50%. In figure 5*b*, an SVPC was generated using the modified algorithm to improve the quality of the lattice. In this case, the beam flowed around both bends without any significant scattering or spreading of the beam. Transmission through the second SVPC approached 97%.

### (b) Multi-mode waveguide bends based on spatially variant band gap materials

More than 10 years ago, a surge of papers emerged on the topic of sharp bends in photonic crystal waveguides, but the vast majority of this work focused on single-mode devices [66–69]. Multi-mode waveguides are of great interest for optical interconnects because different signals can be multiplexed onto different modes, increasing the information capacity [70]. Unfortunately, very little work can be found on the topic of sharp bends in multi-mode waveguides. Some research in multi-mode waveguide bends used TO [70], while others studied photonic crystal waveguides [71–73]. For dense integration, it is desired to form sharp bends without causing inter-modal coupling that would scramble the signals multiplexed onto each mode. Inter-modal coupling remains the primary challenge for multi-mode waveguides. A multi-mode photonic crystal waveguide with a sharp bend was designed in [73] using a topology optimization, but it had limited bandwidth and it is not clear how well the method would extend to waveguides supporting more modes.

A preliminary study on multi-mode photonic crystal waveguides in spatially variant lattices suggests that this is a highly promising technique for realizing sharp bends while also suppressing coupling between the modes. The results from this preliminary study are summarized in figure 6. The devices were simulated using the same finite-difference frequency-domain method identified above. Each row in the figure shows the first three modes in a waveguide propagating around a bend of radius 3.2λ_{0}. The first column shows the first-order mode in each waveguide, the second column shows the second-order mode, and the third column shows the third-order mode. The first row shows a dielectric waveguide in air with a width of 0.5λ_{0} and a core refractive index of 3.4. This is a very high contrast waveguide but showcases dielectric waveguides in their best case for guiding around tight bends. The shape of the modes changes through the bend, indicating that the bend has produced coupling between the modes. The third-order mode is also observed to be leaking out of the waveguide around the bend. The second row shows a uniform photonic crystal with a multi-mode waveguide removed from it which has a bend radius of 3.2λ_{0}. This lattice possesses a band gap so there is no leakage from the waveguides, but the bend induced strong reflections that can be discerned from the field intensity in the plots. In addition, the inter-modal coupling is obvious by looking at the mode profiles exiting the bends which no longer resemble the pure modes. The last row in the figure shows the same photonic crystal, but the lattice was spatially varied to form the bend of radius 3.2λ_{0}. The lattices were generated using deformation control where the critical regions were defined to be the two rows of unit cells on either side of the waveguide. The quality of the lattice outside of these regions is poor, but these regions are not important to performance of the bend. Very clean transmission was observed for all three modes, suggesting the inter-modal coupling was virtually eliminated. While not shown here, similar lattices generated using the basic algorithm without controlling deformations were observed to leak from the waveguide because the deformations weakened the band gap.

## 3. Spatially variant anisotropic metamaterials

Virtually all electric components possess an electromagnetic near-field arising from its reactance. The near-fields are not propagating waves and they remain attached to the component. When two components are placed in close proximity, their near-fields overlap causing the two components to become electromagnetically coupled. This can alter their behaviour and produce crosstalk between the components that degrades signal integrity [74–76]. When a component is embedded in an anisotropic medium, the near-field changes its shape in a manner that prefers to develop in the directions with the highest constitutive values [76]. Positive uniaxial anisotropy is the most straightforward to use here because the near-field prefers only a single direction. When the orientation of positive uniaxial anisotropy is spatially varied in the vicinity of a component, the near-field takes on the shape of the spatial variance [76]. By sculpting the near-field around components, EMC can be improved [76]. Figure 7 illustrates the concept by reshaping the near-field around a thin wire. In figure 7*a*, a thin wire is embedded in an ordinary dielectric medium with relative permittivity set to *ε*_{r}=7.0. The electric field amplitude is symmetric around the wire as expected. In figure 7*b*, the same wire is embedded in a uniform anisotropic medium with *ε*_{xx}=2.0 and *ε*_{yy}=12.0. In this case, the near-field is observed to prefer the *y*-direction because *ε*_{yy}>*ε*_{xx}. In figure 7*c*, the orientation of the anisotropy is spatially varied in the horizontal direction. The near-field is observed to take on the shape of the spatial variance.

The degree to which the near-fields can be sculpted depends on the strength of the anisotropy [76], so it is important to understand how to maximize it. Uniaxial materials exhibit an ordinary permittivity *ε*_{o} in two directions and an extraordinary permittivity *ε*_{e} in the third. The strength of the anisotropy *Δε* can be defined as *Δε*=*ε*_{e}−*ε*_{o}. Positive uniaxial materials have *Δε*>0, while negative uniaxial materials have *Δε*<0. Nature has not provided us with many well-behaved materials that are strongly anisotropic, but this can be achieved using metamaterials. The simplest structure that is positive uniaxial is composed of an array of dielectric rods with high permittivity residing in a background with low permittivity [77]. These metamaterials have very low loss because they do not contain resonant metallic elements. They are also extremely broadband, operating from DC up to a cut-off frequency where the lattice becomes resonant. Using these metamaterials, near-fields can be sculpted on an extremely sub-wavelength scale. Even static fields can be sculpted.

The Weiner bounds given in equation (3.1) set a theoretical limit for the minimum and maximum values of effective permittivity that can achieved [78],
3.1The parameter *f* is the volumetric fill fraction of *ε*_{2} in the metamaterial. An upper limit on the strength of the anisotropy Δ*ε* as a function of fill fraction can be derived from the difference between these limiting values. This is given by
3.2The numerator in equation (3.2) implies that higher contrast lattices will provide stronger anisotropy. Given *ε*_{2}>*ε*_{1}, the fill fraction that maximizes the anisotropy is given by
3.3The higher the contrast between *ε*_{1} and *ε*_{2}, the higher the fill fraction will need to be to maximize the anisotropy. At the optimum fill fraction *f*_{best}, equation (3.2) reduces to equation (3.4) to give the maximum possible anisotropy ,
3.4In practice, usually only 60–70% of this maximum value is realized because the simplest structures do not fully exploit the underlying physics producing the anisotropy.

Figure 8 compares Δ*ε* calculated from the Weiner bounds with the value simulated from a square array of dielectric cylinders in air as a function of fill fraction. The effective tensor quantities were calculated using the plane wave expansion method described in [77]. The simulated structure reaches its maximum Δ*ε* at a lower fill fraction and also achieves a lower Δ*ε* than the theoretical maximum. More complicated structures can push the anisotropy closer to the fundamental limit, but they are likely to be more difficult and expensive to manufacture. The kink observed in the simulated curve between 70% and 80% fill fraction arises because the diameter of the rods exceeds the lattice spacing at this point, causing fill fraction to scale differently.

In prior work, an all-dielectric SVAM was demonstrated to isolate a microstrip transmission line from a metal ball placed in close proximity [76]. The physical mechanism behind the demonstration is illustrated in figure 9*a*. A microstrip transmission line was connected to a vector network analyser that measured the scattering parameters through the line. When a metal ball was placed in close proximity to the microstrip, the scattering parameters fluctuated in excess of 6 dB. To mitigate this sensitivity, an all-dielectric SVAM was designed to tilt the near-field away from the ball so that the ball's presence would not be felt by the microstrip. It was composed of a square array of holes formed in ABS plastic where the orientation of the holes tilted away from the position of the metal ball. TO was circumvented for this design by constructing the input maps to the synthesis algorithm directly. The holes were aligned vertically near the ends of the device and tilted to 60° in proximity to the location of the metal ball. The tilt followed a simple Gaussian profile to conform around the ball. The holes were backfilled with titanium dioxide powder to provide a higher permittivity than the ABS. A photograph of the SVAM used in this experiment is provided in figure 9*b*. After the SVAM was in place, fluctuations of the scattering parameters were reduced to well below 0.5 dB. This is the first known demonstration of an SVAM to improve EMC.

## 4. Conclusion

To make periodic structures that are macroscopically inhomogeneous, it is necessary to spatially vary different geometric attributes of the lattice without also deforming the unit cells and without introducing defects. This is a necessary step in any design process based on spatial transforms. The technique presented in appendix A can spatially vary, or functionally grade, all of the attributes of a periodic structure independently and simultaneously while still producing a lattice that is smooth, continuous and free of defects. Attributes include the structure within the unit cell, orientation of the unit cells, lattice spacing, fill fraction, material composition and even the symmetry of the lattice. In addition to complementing TO, the technique itself enables entirely new device concepts that cannot be designed using TO. First, spatially variant self-collimating photonic crystals were demonstrated to flow light around a very tight bend using ordinary materials with low refractive index. Second, a multi-mode photonic crystal waveguide flowed light around a tight bend without coupling power between the modes, thus maintaining the modal configuration around the bend. Third, SVAMs were shown to have the ability to sculpt the near-field around electric components. This was demonstrated to improve EMC by isolating a microstrip transmission line from a metal ball placed in close proximity.

There are many directions for future work in this area. First, more work is needed to implement the algorithm described in appendix A on high-performance computing platforms [79,80] in order to generate larger and more complex lattices. More work is needed in connecting the tool with TO in order to build devices with arbitrary shapes. Second, the algorithm holds great promise for designing aplanatic metasurfaces. The curvature of these devices can be used to reduce aberrations [81], relax the properties required of the unit cells [82] or simply just to conform to a curved surface such as a radome [83–85]. Third, spatially variant self-collimation is in its infancy. More research is needed to optimize performance, to study the properties of a beam as it propagates through spatially variant discontinuities and to insert the technology into applications such as integrated optics [86] and optical interconnects [87]. Research is needed to study more rigorously the inter-modal coupling in multi-mode SVPC waveguides. This work should include waveguides supporting more modes and waveguides formed from hexagonal lattices supporting wider band gaps. Last, more research is needed to study the simultaneous decoupling of multiple components in close proximity using SVAMs. This will find applications in signal integrity [88–90], EMC [91,92], multiple-input and multiple-output antennas [93–95], phased array antennas [96] and more.

## Authors' contributions

R.C.R. was the primary author of the manuscript, and developed and tested the spatially variant synthesis tool including deformation control. J.J.P. performed the design and simulation of devices in this manuscript; J.L.D. fabricated and characterized the infrared beam-bending SVPCs; S.M.K. developed the DLW process and modifications for smaller feature sizes.

## Competing interests

We declare we have no competing interests.

## Funding

This research was funded in part by DARPA YFA N66001-11-1-415 and NSF CAREER Award DMR/CHE-0748712.

## Acknowledgements

The authors acknowledge the facilities and expertise in three-dimensional printing offered by the W. M. Keck Center for 3D Innovation at the University of Texas at El Paso.

## Appendix A. Improved algorithm for generating spatially variant lattices

A periodic structure is formed any time the material properties are changed as a function of position and the changes repeat themselves in some manner. For electromagnetics, the material properties include the permeability, permittivity and/or conductivity. The following discussion treats only the permittivity, but the technique is easily extended to spatially vary any combination of the constitutive parameters. The standard algorithm for spatially varying periodic structures was first described in [44]. This appendix briefly summarizes this algorithm and then goes on to describe several improvements that increase its efficiency and generate better quality lattices. The basic algorithm is described in two steps. First, the procedure for spatially varying planar gratings is described. Second, this procedure is then generalized to spatially vary arbitrary lattices.

**(a) Basic algorithm to generate spatially variant planar gratings**

The permittivity function *ε*_{a}(** r**) for a planar grating is characterized by its complex amplitude

*a*and its grating vector

**according to equation (A 1), where**

*K***is the position vector, A 1The grating vector conveys two pieces of information at the same time, the direction of the grating and its period. The grating vector is defined to be perpendicular to the grating planes and its magnitude is defined to be 2**

*r**π*divided by the spacing between adjacent planes.

If the period and/or direction of the grating is to be varied as a function of position, then the grating vector ** K**(

**) becomes a function of position and describes the spatial variance of these two parameters at the same time.**

*r***(**

*K***) is called the**

*r**K*-function of the spatially variant grating. The problem is that equation (A 1) fails to construct the grating correctly when the grating vector is a function of position [44]. In order to do this correctly, an intermediate function called the grating phase

*Φ*(

**) is introduced. The grating phase is related to the**

*r**K*-function through the following equation: A 2To solve numerically for the grating phase, equation (A 2) can be cast into matrix form using a numerical technique such as the finite-difference method [44,97]. In this equation, the grating phase is a single scalar quantity that attempts to control all three components of the

*K*-function. There are not enough degrees of freedom to do this perfectly so the solution is a best fit and is calculated using a method such as least squares [44,97]. Given the grating phase, the permittivity is calculated directly from it using equation (A 3), A 3Equations (A 2) and (A 3) are used in place of equation (A 1) to produce a smooth analogue grating. However, the grating does not perfectly correspond to the original

*K*-function. It is a best fit where the errors in the period and orientation of the grating are minimized with equal preference throughout the entire lattice.

In most cases, it is desired to generate binary gratings so that they can be made from a single material. The binary grating *ε*_{b}(** r**) is calculated from the analogue grating using a threshold function

*γ*(

**) according to equation (A 4), A 4All points in the analogue grating with values less than or equal to the threshold are assigned a new value of**

*r**ε*

_{1}. All other points are assigned a value of

*ε*

_{2}. It is possible to generalize equation (A 4) to use multiple threshold values in order to generate lattices containing more than two types of materials. Metallo-dielectric structures containing air, dielectric and metals can be generated this way. Making the threshold a function of position allows the fill fraction to be spatially varied throughout the lattice in addition to the direction and period of the grating. The fill fraction can be spatially varied exactly, so a best fit is not necessary for this parameter. The specific values for the threshold to realize a desired fill fraction can be calculated either numerically or analytically, but they are derived from the geometry of the unit cells. Given the desired fill fraction

*f*(

**) for**

*r**ε*

_{2}, the threshold value for a sinusoidal grating can be estimated as A 5Figure 10 illustrates this process and demonstrates its versatility. In the first column of this figure, the orientation of a binary grating is spatially varied so that it is orientated at 45° within a heart-shaped region and oriented at 0° outside of the heart. The second column shows a binary grating where the fill fraction is spatially varied in the shape of the letter ‘E’. In the third column, the period of the grating is chirped diagonally from the top-left to the bottom-right. All three patterns are applied at the same time in the last column of this figure to demonstrate that all of the attributes of a lattice can be spatially varied at the same time and in unique patterns. In each of these cases, the overall lattice is smooth, continuous and free of defects despite the spatial variance.

A block diagram of the algorithm to spatially vary a planar grating is summarized in figure 11. The input to the algorithm is maps of the desired grating orientation *θ*(** r**), fill fraction

*f*(

**) and the grating period**

*r**Λ*(

**). From these, the target**

*r**K*-function

*K*_{T}(

**) is calculated according to A 6Next, the grating phase is calculated from the**

*r**K*-function by solving equation (A 2). The analogue grating is then calculated from the grating phase using equation (A 3). The last step is to calculate the binary grating by applying the threshold function in equation (A 4). The threshold function

*γ*(

**) is determined from**

*r**f*(

**). The precise relationship depends on the geometry of the unit cell. The output of the algorithm is a spatially variant grating that is the best compromise between grating orientation and grating period. No compromises are needed to spatially vary the fill fraction.**

*r***(b) Modified algorithm to control deformations**

The basic algorithm described above does an excellent job of incorporating all of the spatial variance while minimizing unintentional deformations to the grating. However, some deformations may still remain because the geometry was calculated from a best fit. When there exists additional design freedom, the quality of the lattice can be greatly improved. For example, the anisotropic metamaterials discussed in §3 are highly robust to fluctuations in the lattice spacing. This freedom can be exploited to improve the orientation of the unit cells at the cost of the lattice spacing. It is also possible to improve the lattice spacing at the cost of the orientation. Figure 12 compares a grating generated by the basic algorithm to two other gratings generated with modified algorithms. Figure 12*a* shows the grating generated by the basic algorithm. It possesses the best compromise between the orientation and the period of the grating. Figure 12*b* shows the grating generated with an algorithm modified to enforce only the orientation. In this case, the orientation of the grating is nearly perfect throughout, but the period varies considerably. Figure 12*c* shows the grating generated with an algorithm modified to enforce only the period. The heart pattern is barely discernible and the orientation of the grating is slightly different from 45° within the heart region, but the period is nearly perfect throughout.

The manner in which the algorithm is modified to control deformations is summarized in figure 13. The modified algorithm begins by constructing the target *K*-function *K*_{T}(** r**), which represents the perfect grating. This becomes the first

*K*-function

*K*_{i}(

**) for the iterative algorithm. The grating phase**

*r**Φ*

_{i}(

**) for the first iteration is calculated from**

*r*

*K*_{i}(

**) by solving equation (A 2) as previously discussed. This calculation comes from a best fit so**

*r**Φ*

_{i}(

**) does not correspond to**

*r*

*K*_{i}(

**) exactly. The actual**

*r**K*-function that results is calculated from

*Φ*

_{i}(

**) as**

*r*

*K*_{a,i}(

**)=∇**

*r**Φ*

_{i}(

**). If**

*r*

*K*_{a,i}(

**) were used as the input for a second iteration, the grating phase from the second iteration would be the same as the first. Instead, the**

*r**K*-function for the next iteration

*K*_{i+1}(

**) is constructed from some combination of**

*r*

*K*_{T}(

**) and**

*r*

*K*_{a,i}(

**), depending on what is being enforced. The three simplest approaches for combining these terms are described here. First, to enforce the orientation of the grating at the cost of the period, the new**

*r**K*-function is calculated according to . This calculation retains the magnitude information from the previous iteration, but replaces the angle with the angle component of the target

*K*-function. Second, to enforce the grating period at the cost of the orientation, the new

*K*-function is calculated according to . This retains the angle information from the previous iteration, but replaces the magnitude with the magnitude component of the target

*K*-function. Third, to improve the quality of the grating in critical regions at the cost of the non-critical regions,

*K*_{i+1}(

**) is set equal to**

*r*

*K*_{T}(

**) in the critical regions and set equal to**

*r*

*K*_{a,i}(

**) otherwise. After the new**

*r**K*-function is constructed, a new grating phase function is calculated and the process repeats. Iteration terminates when the change in the grating phase from one iteration to the next is sufficiently small. This change can be calculated according to equation (A 7), where

*S*is the volume or area of the lattice, A 7In this work, iteration continued until Δ

*Φ*

_{i+1}<2×10

^{−5}. This happened within 20–40 iterations, but the number can change depending on the size of the lattice and the complexity of the spatial variance. Hybrids of these three approaches are possible along with using weighted combinations of the target and actual

*K*-functions.

**(c) Spatially varying arbitrary lattices**

The permittivity of a periodic structure is a periodic function, so it can be expanded into a complex Fourier series as in equation (A 8),
A 8
A 9
A 10Comparing equation (A 8) with equation (A 1), the term inside the summation has the form of a planar grating so it is concluded that a lattice can be decomposed into a set of planar gratings of the same physical size. The complex amplitudes of the planar gratings are calculated by Fourier transforming the unit cell. Numerically, this can be done using a fast Fourier transform as indicated by equation (A 9). The grating vectors for each planar grating are calculated analytically from the reciprocal lattice vectors *T*_{1}, *T*_{2} and *T*_{3} [98] using equation (A 10). The terms *p*, *q* and *r* are integers identifying the different planar gratings in the expansion. The algorithm runs faster when fewer planar gratings are retained, but the unit cells can be distorted depending on the severity of the truncation.

A block diagram of the algorithm to spatially vary arbitrary lattices is provided in figure 14. It begins with the same input information, but now also includes a description of the unit cell. Greyscale unit cells are better than binary unit cells when the fill fraction is to be adjusted at the end. The first task in the algorithm is to decompose the unit cell into its component planar gratings. This involves calculating the Fourier coefficients *a*_{pqr} and the corresponding grating vectors *K*_{pqr}(** r**). From here, it is most efficient to minimize the number of planar gratings that must be processed. One thing that can be done is to ignore all planar gratings with Fourier coefficients

*a*

_{pqr}that are less than some small threshold. This usually achieves a 60–90% reduction in the number of planar gratings without significantly distorting the unit cells. A second thing is to sort the remaining planar gratings into sets containing all parallel grating vectors. For each set, only a single solution is needed from the algorithm summarized in figure 13. The grating phase for all other gratings in the set can be calculated immediately from this single solution using equation (A 11), A 11This second approach usually achieves a 70% reduction for two-dimensional lattices, and a 60% reduction for three-dimensional lattices. Using both together easily achieves a 90–99% overall reduction. Visualizations of these truncation schemes are provided at each step in figure 14. Whatever planar gratings remain in the expansion are each spatially varied over the volume of the lattice as described previously. Deformation control is applied to each planar grating separately. The overall analogue lattice is computed by summing all of the planar gratings according to A 12The ‘

*Re*’ operation in this equation extracts only the real part from the sum in order to remove small residual imaginary components that originate because the solution is obtained numerically.

The analogue lattice can be analysed to quantify any deformations to the lattice that may still exist. In many circumstances, it is possible to electromagnetically compensate for these stubborn deformations by spatially varying another property of the lattice. For example, when the lattice spacing is compressed, the fill fraction can be raised to increase the electrical size of the unit cells. When the lattice spacing is stretched, the fill fraction can be lowered to reduce the electrical size of the unit cells. The manner in which the fill fraction is adjusted to compensate must come from a rigorous electromagnetic simulation of the deformed unit cell. Given this map of desired fill fraction, the spatially variant threshold function is calculated. The final step in the process is to apply the threshold function to spatially vary the fill fraction. This compensation step is subtle, but illustrated in the bottom part of figure 14.

Figure 15 compares spatially variant lattices where there exists some critical region where the geometry of the lattice is most important. This critical region is highlighted in red and this may be the desired path for a beam. Figure 15*a* shows a lattice generated using the basic algorithm summarized in figure 11. Large fluctuations in the lattice spacing are observed in the critical regions. Figure 15*b*,*c* shows the same lattice generated using 10 and 40 iterations of the modified algorithm summarized in figure 13. The lattice inside the critical region is greatly improved at the cost of the lattice outside of the critical region. This is acceptable because no electromagnetic waves should be present outside of the critical region. The algorithm was observed to converge for this lattice at around 40 iterations.

**(d) Spatially varying arrays of discontinuous metallic elements**

It is inefficient to decompose thin metallic elements into a set of planar gratings because many hundreds or thousands of planar gratings would be needed in the expansion in order to represent the elements accurately. Instead, it is possible to generate these lattices using just two or three planar gratings. For two-dimensional lattices and metasurfaces, only two planar gratings are needed and these are spatially varied within the plane of the device. Figure 16*a* shows the superposition of two spatially variant planar gratings along with the outline of a device generated from this map. Metallic elements are placed at the intersections of the two planar gratings in an orientation defined by the planar gratings at the intersections. An example lattice generated from these data is shown in figure 16*b*. The specific geometry of the metallic elements at each point comes from the effective properties that they must realize at these points. This is essentially step 3 in figure 1. The position and orientation of the metallic elements comes from the spatially variant planar gratings.

**(e) Periodic structures on curved surfaces**

It is sometimes desired to put a periodic structure such as a metasurface onto a curved substrate. Such aplanatic designs are difficult to generate without deforming the lattice due to the curvature. Even when it is not desired to functionally grade the metasurface, the lattice must still be spatially varied in order to keep it uniform over the curved geometry. An example of this process is given in figure 17 where a hexagonal array of metal rings was put onto a curved surface uniformly. The algorithm for doing this is the same as discussed above except that the directions in which the derivatives are calculated must be modified to conform to the curvature. When this modification is incorporated, equation (A 2) becomes A 13and A 14

A simple way to solve equation (A 13) is with the finite-difference method. One way to do this is to mesh the surface like that shown in figure 17*a*. Here, the *x*- and *y*-coordinates of the vertices are distributed evenly and it is only the *z*-coordinate that is changed to describe the surface height variations. More sophisticated meshing is certainly possible. The new directions that the derivatives are calculated are indicated by the red and green arrows superimposed onto the mesh. What used to be the *x*-direction in equation (A 2) is now the direction in equation (A 13), and what used to be the *y*-direction is now the direction. These new directions are chosen to conform to the curved surface and are always tangential to it. The derivatives in equation (A 4) can be approximated with central finite-differences according to equation (A 15), where *r*_{p,q} is the physical position of the grid vertex identified by the array indices *p* and *q*,
A 15

The two *K*-functions needed to construct a hexagonal lattice are now defined to be in the and directions. If the lattice spacing is *Λ*, these are calculated according to
A 16

The remainder of the algorithm is the same including deformation control, placing discrete metal elements, and spatially varying arbitrary lattices. Figure 17*c*,*d* shows two views of the two planar gratings generated from the *K*-functions in equation (A 16). Figure 17*e* shows these two planar gratings superimposed to form a hexagonal lattice that is evenly distributed over the curved surface. The last step for a metasurface is to place metal elements at the intersections of the gratings as shown in figure 17*b*. In this case, the elements are metal rings and they are oriented to be parallel to the surface.

## Footnotes

One contribution of 14 to a Theo Murphy meeting issue ‘Spatial transformations: from fundamentals to applications’.

- Accepted February 12, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.