## Abstract

It is often desired to functionally grade and/or spatially vary a periodic structure like a photonic crystal or metamaterial, yet no general method for doing this has been offered in the literature. A straightforward procedure is described here that allows many properties of the lattice to be spatially varied at the same time while producing a final lattice that is still smooth and continuous. Properties include unit cell orientation, lattice spacing, fill fraction, and more. This adds many degrees of freedom to a design such as spatially varying the orientation to exploit directional phenomena. The method is not a coordinate transformation technique so it can more easily produce complicated and arbitrary spatial variance. To demonstrate, the algorithm is used to synthesize a spatially variant self-collimating photonic crystal to flow a Gaussian beam around a 90° bend. The performance of the structure was confirmed through simulation and it showed virtually no scattering around the bend that would have arisen if the lattice had defects or discontinuities.

©2012 Optical Society of America

## 1. Introduction

It is ultimately desired to control everything about the electromagnetic field throughout some volume. This includes the flow of energy, amplitude profile of the field, phase of the waves, polarization, coupling between modes, and potentially more. Many degrees of freedom are needed in order to control this many properties. Periodic structures like photonic crystals [1,2] and metamaterials [3] have been shown to provide extraordinary control over the electromagnetic field, offering many new and novel functionality that cannot be obtained using naturally occurring materials. Further, light cannot be controlled using homogeneous structures. Devices must be composed of interfaces, graded index, or other inhomogeneous features. To meet this demand, a general purpose synthesis procedure was developed that is capable of simultaneously spatially varying the orientation, lattice spacing, fill fraction, and other properties of a periodic structure throughout its volume in a way that leaves the overall lattice smooth and continuous. Avoiding discontinuities in the lattice is important because these produce scattering and field concentrations that can degrade the performance of the device. The ability to spatially vary the orientation of the unit cells throughout a lattice enables directional phenomena like band gaps, anisotropy, and dispersion to be fully exploited. Having a tool capable of varying so many degrees of freedom is critical when multiple properties of the field and devices are to be controlled at the same time. The synthesis algorithm described here provides this freedom.

The method outlined in this paper is unique in that it calculates the physical geometry of a spatially variant lattice so that it can be directly fabricated. In contrast, transformation optics is a coordinate transformation technique that only calculates the permittivity and permeability tensor functions [4]. It does not calculate the physical geometry of the lattice that is needed to realize the material tensors. This information must come from a subsequent step and the synthesis tool can serve this function. In other work, graded refractive index profiles were realized by spatially varying the fill fraction of a lattice [5–9]. The synthesis tool described here is capable of doing this and much more. In fact, virtually any set of properties of the lattice can be spatially varied at the same time while still rendering the final lattice smooth and continuous.

This paper begins by describing the synthesis algorithm in detail, providing several examples of spatially variant lattices and some notes on efficient implementation. The technique is demonstrated by synthesizing and simulating a spatially variant self-collimating photonic crystal that flows a Gaussian beam around a 90° bend with virtually no scattering introduced by the bend. To the author’s knowledge, this is the first known example of spatially variant lattice orientation to control the flow of energy using self-collimation. This approach could provide an attractive alternative to transformation optics and graded index devices because the final structure requires no anisotropy, permeability, or extreme parameter values and the approach provides many additional degrees of freedom to control the field. A finite-difference time-domain (FDTD) simulation of the lattice showed virtually no scattering of the beam induced by the bend. This was achieved because the synthesized lattice was smooth and continuous. While self-collimation was demonstrated here, the synthesis algorithm is broadly applicable to a wide range of phenomena and applications.

## 2. Synthesis algorithm

The following sections describe the synthesis algorithm in two parts. Section 2.1 discusses the data required by the synthesis algorithm including the baseline unit cell and how it is to be spatially varied throughout a device. Section 2.2 describes the actual synthesis algorithm that processes the input data to generate the final spatially variant lattice.

#### 2.1 Construct input data

Before the synthesis algorithm can be implemented, several functions must be constructed that define the baseline unit cell and how it is to be spatially varied as a function of position. These are inputs to the synthesis procedure and are not generated by the algorithm itself.

#### Design the unit cell and compute its spatial harmonics

The baseline unit cell describes the specific geometry that is to be arrayed and spatially varied to form the final lattice. It must be designed prior to the synthesis procedure and is typically derived from a rigorous electromagnetic optimization to produce whatever property is desired. To keep the description of the algorithm generic and more illustrative, a simple unit cell composed of a centered triangle will be used. In this paper, the unit cell function _{${\epsilon}_{\text{uc}}\left(\overrightarrow{s}\right)$} is initialized to all zeros and then set to 1 anywhere that dielectric is to be placed. This can be extended to incorporate multiple materials. The actual dielectric constants of the materials are incorporated in a later step.

Assuming the unit cell is periodic, it can be decomposed into a complex Fourier series. Each term in the Fourier series is called a spatial harmonic and can be interpreted as a one-dimensional (1D) sinusoidal grating that when summed together reconstruct the original unit cell. Each of these has its own unique period, complex amplitude, and direction. While an infinite number of spatial harmonics is needed to construct the unit cell exactly, the series must be truncated to a finite set of *M* spatial harmonics for practical implementation on a computer. While the choice of which spatial harmonics are retained is arbitrary, typically all of the lowest order harmonics up to some cutoff are selected.

In this equation, _{$\overrightarrow{s}$} is position, *a _{m}* is the complex amplitude of the

*m*

^{th}spatial harmonic, and

_{${\overrightarrow{K}}_{m}$}is the grating vector associated with the

*m*

^{th}spatial harmonic. The complex amplitudes can be calculated using a fast Fourier transform (FFT). The grating vectors associated with the spatial harmonics for orthorhombic symmetries are calculated using Eq. (2) where

*p*,

*q*, and

*r*are integers and Λ

*, Λ*

_{x}*, and Λ*

_{y}*are the unit cell dimensions in the*

_{z}*x*,

*y*, and

*z*directions respectively.

The concept of decomposing a unit cell into its spatial harmonics is illustrated in Fig. 1 for both two-dimensional (2D) and three-dimensional (3D) unit cells. In Fig. 1(a) is the 2D unit cell that will be used as an example throughout this section. The lowest order 5 × 5 spatial harmonics are shown as 1D sinusoidal gratings calculated over the area of the unit cell. When all of the spatial harmonics are superimposed, the sum reconstructs the original unit cell. In Fig. 1(b) is a 3D unit cell along with its lowest order 5 × 5 × 5 spatial harmonics.

As an alternative to above, periodic structures can be described using very few spatial harmonics when the grating vectors are set equal to some linear combination of the reciprocal vectors of the lattice [10]. Using a small number of spatial harmonics limits the geometry of the unit cells, but speeds the synthesis procedure and enables a spatially variant fill fraction to be incorporated more easily. Borrowing from holographic lithography, all 14 Bravais lattice symmetries can be realized using just three grating vectors [11]. Lattices uniform in one direction can be realized with just two grating vectors and lattices uniform in two directions can be realized with just one grating vector. Common formulations for the grating vectors for cubic symmetries are provided in Fig. 2 .

#### Define the spatially variant parameters

To control how the unit cell orientation, lattice spacing, and fill fraction are to be spatially varied throughout the volume of the lattice, separate functions must be constructed that define each of these parameters as a function of position. These are defined prior to the synthesis procedure to make the final device perform whatever task is desired. Examples of these functions and their individual effects on the resulting lattice are shown in Fig. 3 . Figure 3(h) shows the final lattice when all of the parameters are spatially varied at the same time. Describing the spatial variance in this manner is more intuitive and offers the ability to incorporate much more complicated patterns than coordinate transformation techniques. In addition, each function can defined completely independently.

#### 2.2 Synthesize the spatially variant lattice

The synthesis algorithm generates the final lattice by spatially varying each of the spatial harmonics separately and then adding the results. The steps to do this for each spatial harmonic are described below.

#### Step 1: Construct spatially variant K function

A grid is constructed over the entire volume of the lattice and a spatially variant grating vector function (*K* function) is computed throughout this grid as illustrated in Fig. 4
.First, the grating vector _{${\overrightarrow{K}}_{m}$} associated with the *m*^{th} spatial harmonic is computed from Eq. (2). This is then distributed uniformly across the grid to arrive at the uniform *K* function. Second, the tilt of the orientation function is added to this *K* function to arrive at an intermediate *K* function. Third, the magnitude of the *K* function is divided by the lattice spacing function to produce the final spatially variant *K* function. In the end, the spatially variant *K* function describes both the orientation and period of the 1D grating as a function of position throughout the lattice.

#### Step 2: Calculate the grating phase function

Having seen Eq. (1), it is tempting to construct the 1D grating associated with the *m*^{th} spatial harmonic directly from the spatially variant *K* field according to the following equation.

*K*function is spatially varied. Instead, the spatially variant 1D grating can be constructed using an intermediate function called the grating phase. The grating phase

_{${\varphi}_{m}\left(\overrightarrow{s}\right)$}is related to

_{${\overrightarrow{K}}_{m}\left(\overrightarrow{s}\right)$}through a gradient operation.

A number of analytical and numerical methods can be used for solving Eq. (4) to calculate the grating phase from the spatially variant *K* function. In this work, the finite-difference method [12] was used to approximate the derivatives in Eq. (4). This led to a large set of equations that were written in matrix form as

The term **Φ*** _{m}* is a column vector containing the grating phase at each point on the grid. These values are not known at this point, but will be calculated shortly. The terms

**D**

*,*

_{x}**D**

*, and*

_{y}**D**

*are sparse banded matrices that perform derivative operations across the grid using central finite-differences. The terms*

_{z}**K**

*,*

_{x,m}**K**

*, and*

_{y,m}**K**

*are column vectors containing the components of the*

_{z,m}*K*function across the grid.

This matrix formulation contains more equations than it does unknowns. This means there are more attributes of the grating that are trying to be controlled than there are degrees of freedom to do so. For this reason, the matrix equation is solved in the sense of least squares [13] and the resulting 1D grating is seen as a best fit. The solution in the sense of least squares is summarized in Eqs. (6)–(8). The matrix _{${A}^{\prime}$} remains a sparse matrix and Eq. (6) should be solved as a sparse system for best efficiency [14].

In Fig. 5 , direct construction using Eq. (3) is compared to grating phase construction using Eqs. (6)–(8). The 1D grating produced by direct construction is shown in Fig. 5(a). In this case, the lattice spacing and orientation are greatly distorted. The 1D grating constructed using the grating phase approach is shown in Fig. 5(b). It has minimized any potential lattice distortions.

#### Step 3: Compute spatially variant 1D grating

Given the grating phase function _{${\varphi}_{m}\left(\overrightarrow{s}\right)$} and the complex amplitude *a _{m}* of the

*m*

^{th}spatial harmonic, it is straightforward to compute the corresponding spatially variant 1D grating using Eq. (9). While not discussed here, the amplitude terms can be spatially varied to change the unit cell geometry throughout the lattice. The dielectric function

_{${\epsilon}_{m}\left(\overrightarrow{s}\right)$}is in general complex to convey the offset of the grating. While it is not needed at this point, the 1D grating can be visualized by ignoring the imaginary component.

#### Step 4: Construct overall lattice

Given all of the spatially variant 1D gratings, a preliminary spatially variant lattice _{${\epsilon}^{\prime}\left(\overrightarrow{s}\right)$} is constructed from their sum. It is preliminary because the spatially variant fill fraction has not yet been incorporated.

Due to numerical artifacts arising from the FFT and construction of the lattice using Eq. (10), the resulting dielectric distribution may still have a small complex component. This is typically many orders of magnitude less than the real part and it should be neglected.

#### Step 5: Incorporate spatially variant fill fraction

The spatially variant fill fraction of the lattice can be incorporated by applying a threshold technique to the preliminary dielectric distribution according to Eq. (11). The preliminary dielectric distribution was made binary for this reason. This step also produces the final lattice by incorporating the actual dielectric constants *ε*_{1} and *ε*_{2}. It is straightforward to extend this approach to handle more than two materials.

When only a small number of spatial harmonics are used, the preliminary dielectric distribution calculated using Eq. (10) varies in a near sinusoidal manner. In this case, the threshold needed to realize the fill fraction *f* of *ε*_{2} can be estimated according to

The manner in which the fill fraction of the lattice is controlled using a threshold is demonstrated in Fig. 6
where a different lattice with lower contrast is used as the example. As the threshold value *γ* is lowered, the fraction of the lattice filled by *ε*_{2} increases. When many spatial harmonics are used, the resulting dielectric distribution will have high contrast and direct application of Eq. (11) will have little effect. Some control over fill fraction can still be obtained using a grayscale unit cell or by partially blurring the preliminary dielectric function _{${\epsilon}^{\prime}\left(\overrightarrow{s}\right)$} before applying Eq. (11).

#### 2.3 Notes on implementation

The main computational bottleneck in the synthesis procedure is calculating the grating phase using Eq. (6). This is especially true for large grids. Inspection of Fig. 5 shows that the grating phase function varies very slowly compared to the 1D grating it is used to construct. For this reason, it is more efficient to calculate the grating phase on a coarse grid, interpolate that solution onto a fine grid, and then construct the 1D grating on the fine grid. A rule of thumb for the coarse grid is to use at least several points per unit cell in the overall lattice. Given some integer *N*_{coarse} (4 to 10 is typical) and the lattice constant Λ, the resolution of the coarse grid should be

The resolution of the fine grid must be sufficient to resolve the highest order spatial harmonic which will have the shortest period *L*_{min}. Given some integer *N*_{fine}>10, the resolution of the fine grid should be

It is not possible to spatially vary the orientation of unit cells throughout a lattice without also distorting the lattice in some way. This is most severe when the spatial variance is abrupt relative to the size of the unit cell. This can weaken the desired electromagnetic response or shift the response to a different wavelength. In many cases, the properties can be retuned by incorporating a spatially variant fill fraction or by spatially varying some other property to compensate. In Fig. 7(a) , a lattice is illustrated where the spatially variant orientation has stretched the size of the unit cells at the top of the lattice and compressed the size of the unit cells at the bottom. The electrical size of the unit cells at the top of the lattice can be reduced by decreasing the fill fraction to remove dielectric. The electrical size of the unit cells at the bottom of the lattice can be increased by increasing the fill fraction to add dielectric. This method of compensation is illustrated in Fig. 7(b) where the fill fraction is adjusted in a radial pattern to match the distortion in the lattice. The exact amount by which the fill fraction should be adjusted must come from an electromagnetic model, but the pattern in which it is spatially varied can be determined by the spatially variant lattice itself. The total adjustment was exaggerated in this figure for illustration purposes. In practice, it is usually more subtle.

For greatest clarity, only 2D lattices were illustrated and described above. The synthesis technique can easily create fully 3D lattices as well without any modification of the governing equations. In the next section, a 2D spatially variant self-collimating lattice will be presented. The 3D equivalent to this is shown in Fig. 8 . This lattice possesses the ability to fully self-collimate a Gaussian beam and is presented here only to illustrate the concept of spatially varying a 3D lattice.

## 3. Example – spatially variant self-collimating photonic crystal

Some photonic crystals have the intriguing ability to self-collimate a beam of light without requiring nonlinear material interactions [15–18]. The beam propagates without diffraction or “spreading” when the spatial dispersion is engineered such that all the Fourier components of the beam propagate with the same group velocity. In this manner, energy is forced to flow in a specific direction regardless of the shape, direction, or divergence of the applied beam. This directional phenomenon can be exploited in spatially variant lattices to flow energy in useful ways.

The conditions that must exist for self-collimation to occur can be understood through the isofrequency surfaces. An isofrequency surface is like an index ellipsoid [19] and describes the spatial dispersion of the lattice at a specific frequency. The distance from the origin to a point on the surface is proportional to the effective phase refractive index a wave would experience if it was propagating in that direction. In contrast, energy propagates in a direction normal to the isofrequency surface at that point because group velocity is calculated from the gradient of the lattice dispersion in *k*-space.

In periodic structures, the spatial dispersion ∇* _{k}ω* can be engineered so that portions of the isofrequency surfaces are flat over some span of wave vectors

*k*and over some band of frequencies

*ω*. This is the region where self-collimation is supported. The effect can be designed and optimized by tailoring the symmetry and composition of the unit cell.

A simple self-collimating lattice and its isofrequency contours for the second photonic band were calculated using the plane wave expansion method [20]. These are all depicted in Fig. 9
. This unit cell was selected because it can be constructed using just two spatial harmonics with amplitudes *a*_{1} = *a*_{2} = 1.0. The grating vectors associated with the two spatial harmonics were

The fill fraction of the lattice was made 43.7% by setting the threshold value in Eq. (11) to *γ* = 0.1178. Using a material with dielectric constant 2.5, the center frequency of self-collimation was taken from Fig. 9(b) to be Λ/*λ*_{0} ≅ 0.7039. Given this, the lattice constant can be chosen to produce self-collimation at whatever wavelength is desired. For operation at *λ*_{0} = 1550 nm, the lattice constant should be Λ = 1091 nm.

A 2D spatially variant self-collimating photonic crystal lattice was synthesized from the unit cell described above to flow energy around a 90° bend. This was formed by spatially varying the orientation of the unit cells in an azimuthal pattern. Using self-collimation to control the flow of a beam instead of transformation optics avoids the need for anisotropy, permeability, and extreme parameter values. It further provides more degrees of freedom than analogous devices based on graded index. The lattice size was made 40Λ × 40Λ with a grid resolution of Δ = 0.03*λ*_{0}. The finite-difference time-domain (FDTD) method was used to simulate the lattice using the TM* _{z}* mode [21] and the results are shown in Fig. 10
. A uniaxial perfectly matched layer (UPML) [22] was incorporated into the outer 40 grid points to absorb outgoing waves. The overall grid size was 1288 × 1288 and the simulation was stopped after 3,000 iterations. The lattice was illuminated from the top with a Gaussian beam of width 4

*λ*

_{0}. It is important to reinforce that the beam was turned using spatially variant self-collimation and not waveguiding, band gaps, anisotropy, graded fill fraction, or graded refractive index. That leaves all of these mechanisms as additional degrees of freedom to control other aspects of the field or device.

The simulation showed that the beam remained self-collimated while inside the lattice and that it turned the 90° bend with virtually no scattering induced by the bend. Scattering was observed at the edges of the photonic crystal due to the mismatch between air and the lattice. In addition, the output beam was split into two dominant modes. Self-collimation only controls the flow of energy so phase remained uncontrolled in this design. It was surmised that phase was induced asymmetrically due to the bend which produced the asymmetry in the exiting field. Future work will entail spatially varying additional parameters like lattice constant and fill fraction to further control the phase of the beam through the device.

## 4. Conclusion

This paper described a general purpose algorithm to synthesize spatially variant lattices such as photonic crystals or metamaterials. The procedure is not based on coordinate transformation so it can more easily produce complicated patterns as well as spatially vary multiple properties of the lattice simultaneously. This offers much greater design freedom than other approaches because so many attributes of the lattice can be controlled at the same time. The ability to spatially vary the unit cell orientation is useful for exploiting phenomena with directional dependencies. The technique produces a final lattice that is smooth and continuous. This is important because scattering and field concentrations at discontinuities can degrade device performance. In addition, the lattice spacing remains strikingly uniform when the unit cell orientation is spatially varied. This is important for maintaining consistent properties throughout the lattice.

The procedure was demonstrated by synthesizing a spatially variant self-collimating photonic crystal designed to flow a Gaussian beam around a 90° bend. The performance of the device was confirmed using FDTD and the beam showed virtually no scattering around the bend. While self-collimation was demonstrated here, the synthesis tool is broadly applicable to many phenomena and applications. For example, diffraction gratings can be spatially varied to form novel beam shaping and spectral filtering devices or to operate on curved surfaces.

## Acknowledgments

This work was supported in part by a DARPA Young Faculty Award (grant no. N66001-11-1-415).

## References and links

**1. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic Crystals, Molding the Flow of Light* (Princeton University Press, 1995).

**2. **H. Benistry, V. Berger, J.-M. Gerard, D. Maystre, and A. Tchelnokov, *Photonic Crystals, Towards Nanoscale Photonic Devices* (Springer, 2005).

**3. **S. A. Ramakrishna and T. M. Grzegorczyk, *Physics and Applications of Negative Refractive Index Metamaterials*, (CRC Press, 2009).

**4. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science **312**(5781), 1780–1782 (2006). [CrossRef] [PubMed]

**5. **D. H. Spadoti, L. H. Gabrielli, C. B. Poitras, and M. Lipson, “Focusing light in a curved-space,” Opt. Express **18**(3), 3181–3186 (2010). [CrossRef] [PubMed]

**6. **E. G. Johnson, M. K. Poutous, Z. A. Roth, P. Srinivasan, A. J. Pung, and Y. O. Yilmaz, “Advanced fabrication methods for 3D meta-optics,” Proc. SPIE **7927**, 792706, 792706-7 (2011). [CrossRef]

**7. **Z. A. Roth, P. Srinivasan, M. K. Poutous, A. Pung, R. C. Rumpf, and E. G. Johnson, “Azimuthally varying guided mode resonance filters,” Micromachines **3**(1), 180–193 (2012). [CrossRef]

**8. **E. Akmansoy, E. Centeno, K. Vynck, D. Cassagne, and J.-M. Lourtioz, “Graded photonic crystals curve the flow of light: An experimental demonstration by the mirage effect,” Appl. Phys. Lett. **92**(13), 133501 (2008). [CrossRef]

**9. **B. Vasić, G. Isić, R. Gajić, and K. Hingerl, “Controlling electromagnetic fields with graded photonic crystals in metamaterial regime,” Opt. Express **18**(19), 20321–20333 (2010). [CrossRef] [PubMed]

**10. **R. C. Rumpf, “Design and optimization of nano-optical elements by coupling fabrication to optical behavior,” Ph.D. Dissertation, University of Central Florida (2006), pp. 171–183.

**11. **L. Z. Cai, X. L. Yang, and Y. R. Wang, “All fourteen Bravais lattices can be formed by interference of four noncoplanar beams,” Opt. Lett. **27**(11), 900–902 (2002). [CrossRef] [PubMed]

**12. **S. C. Chapra and R. P. Canale, *Numerical Methods for Engineers with Software and Programming Applications*, 4th Ed., 820–856 (McGraw-Hill, 2002).

**13. **B. Noble and J. W. Daniel, *Applied Linear Algebra*, 3rd ed. (Prentice Hall, 1988), pp. 66–73.

**14. **Y. Sadd, *Iterative Methods for Sparse Linear Systems*, 2nd ed. (Yousef Sadd, 2000).

**15. **H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Tamamura, T. Sato, and S. Kawakami, “Self-collimating phenomena in photonic crystals,” Appl. Phys. Lett. **74**(9), 1212–1214 (1999). [CrossRef]

**16. **R. Iliew, C. Etrich, and F. Lederer, “Self-collimation of light in three-dimensional photonic crystals,” Opt. Express **13**(18), 7076–7085 (2005). [CrossRef] [PubMed]

**17. **J. Shin and S. Fan, “Conditions for self-collimation in three-dimensional photonic crystals,” Opt. Lett. **30**(18), 2397–2399 (2005). [CrossRef] [PubMed]

**18. **Z. Lu, S. Shi, J. A. Murakowski, G. J. Schneider, C. A. Schuetz, and D. W. Prather, “Experimental demonstration of self-collimation inside a three-dimensional photonic crystal,” Phys. Rev. Lett. **96**(17), 173902 (2006). [CrossRef] [PubMed]

**19. **M. Born and E. Wolf, *Principles of Optics*, 6th Ed., 673–678 (Cambridge University Press, 1980).

**20. **R. C. Rumpf, “Design and optimization of nano-optical elements by coupling fabrication to optical behavior,” Ph.D. Dissertation, University of Central Florida (2006), pp. 109–124.

**21. **A. Taflove and S. C. Hagness, *Computational Electrodynamics, the Finite-Difference Time-Domain Method*, 3rd ed. (Artech House, 2005).

**22. **Z. S. Sacks, D. M. Kingsland, R. Lee, and J.-F. Lee, “A Perfectly Matched Anisotropic Absorber for Use as an Absorbing Boundary Condition,” IEEE Trans. Antenn. Propag. **43**(12), 1460–1463 (1995). [CrossRef]