Modeling and optimization of photonic crystal devices based on transformation optics method

In this paper, we propose a method for designing Photonic Crystal (PhC) devices that consist of dielectric rods with varying size. In the proposed design method, PhC devices are modeled with the Transformation Optics (TO) approach, and then they are optimized using the gradient method. By applying the TO technique, the original device model is transformed into an equivalent model that consists of uniform and fixed-sized rods, with parameterized permittivity and permeability distributions. Therefore, mesh refinement around small rods can be avoided, and PhC devices can be simulated more efficiently. In addition, gradient of the optimization object function is calculated with the Adjoint-Variable Method (AVM), which is very efficient for optimizing devices subject to multiple design variables. The proposed method opens up a new avenue to design and optimize a variety of photonic devices for optical computing and information processing. © 2014 Optical Society of America OCIS codes: (130.3120) Integrated optics devices; (130.5296) Photonic crystal waveguides; (350.4238) Nanophotonics and photonic crystals. References and links 1. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Ed. (Princeton University, 2008). 2. I. A. Sukhoivanov, and I. V. Guryev, Photonic Crystals: Physics and Practical Modeling (Springer-Verlag, Berlin, 2009). 3. A. Adibi, Y. Xu, R. K. Lee, A. Yariv, and A. Scherer. “Guiding mechanisms in dielectric-core photonic-crystal optical waveguides,” Phys. Rev. B 64, 033308 (2001). 4. P. Bienstman, S. Assefa, S. G. Johnson, J. D. Joannopoulos, G. S. Petrich, and L. A. Kolodziejski, “Taper structures for coupling into photonic crystal slab waveguides,” J. Opt. Soc. Am. B 20, 1817–1821 (2003). 5. K. Dossou, L. C. Botten, C. M. Sterke, R. C. McPhedran, A. A. Asatryan, S. Chen, and J. Brnovic, “Efficient couplers for photonic crystal waveguides,” Opt. Commun. 265, 207–219 (2006). 6. A. Mekis and J. D. Joannopoulos, “Tapered couplers for efficient interfacing between dielectric and photonic crystal waveguides,” J. Lightw. Technol., 19(6), 861–865 (2001). 7. O. Ozgun, M. Kuzuoglu, “Software metamaterials: Transformation media based multiscale techniques for computational electromagnetics”, J. Comput. Phys. 236, 203–219 (2013). #202078 $15.00 USD Received 27 Nov 2013; revised 6 Jan 2014; accepted 9 Jan 2014; published 30 Jan 2014 (C) 2014 OSA 10 February 2014 | Vol. 22, No. 3 | DOI:10.1364/OE.22.002725 | OPTICS EXPRESS 2725 8. A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equation,” J. Mod. Opt. 43(4), 773–793 (1996). 9. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). 10. Y. M. Liu and X. Zhang, “Recent advances in transformation optics,” Nanoscale 4(17), 5277–5292 (2012). 11. L. H. Gabrielli, D. Liu, S. G. Johnson, and M. Lipson, “On-chip transformation optics for multimode waveguide bends,” Nat. Commun. 3:1217 doi: 10.1038/ncomms2232 (2012). 12. Q. Wu, J. P. Turpin, and D. H. Werner, “Integrated photonic systems based on transformation optics enabled gradient index devices,” Light: Science & Applications (2012) 1, e38; doi:10.1038/lsa.2012.38. 13. P. Sanchis, P. Villalba, F. Cuesta, A. Håkansson, A. Griol, J. V. Galán, A. Brimont, and J. Marti, “Highly efficient crossing structure for silicon-on-insulator waveguides,” Opt. Lett. 34(18), 2760–2762 (2009). 14. Y. Zhang, S. Yang, A. E. J. Lim, G. Q. Lo, C. Galland, T. Baehr-Jones, and M. Hochberg, “A compact and low loss Y-junction for submicron silicon waveguide,” Opt. Express 21(1), 1310–1316 (2013). 15. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design.” Opt. Express 21(18), 21693–21701 (2013). 16. N. K. Georgieva, S. Glavic, M. H. Bakr, and J. W. Bandler, “Feasible Adjoint Sensitivity Technique for EM Design Optimization,” IEEE Tans. Microwave Theory Tech. 50(12), 2751–2758 (2002). 17. G. Veronis, R. W. Dutton, and S. Fan, “Method for sensitivity analysis of photonic crystal devices,” Opt. Lett. 29(19), 2288–2290 (2004). 18. J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: a high-bandwidth low-loss T-junction waveguide,” J. Opt. Soc. Am. B 22(6), 1191–1198 (2005). 19. COMSOL Multiphysics, http://www.comsol.com/ 20. D. Marcuse, Theory of Dielectric Optical Waveguides (Academic Press, 1974). 21. S. G. Johnson and J. D. Joanopoulos, “Block-iterative frequency-domain methods for Maxwell’s equation in a planewave basis,” Opt. Express 8(3), 173–190 (2001). 22. M. Schmidt, http://www.di.ens.fr/ ̃Emschmidt/Software/minConf.html 23. M. Minkov and V. Savona, “Effect of hole-shape irregularities on photonic crystal waveguides,” Opt. Lett. 37(15), 3108–3110 (2012). 24. V. Savona, “Electromagnetic modes of a disordered photonic crystal,” Phys. Rev. B 83, 085301 (2011).


Introduction
Photonic Crystals (PhCs) are novel artificial materials that allow us to manipulate electromagnetic waves in an unprecedented manner [1,2].A fundamental property of PhCs is the Photonic Band Gap (PBG) [1,2], which prohibits light propagation in PhCs at certain frequency ranges and can be exploited to control light propagation.By introducing defects in PhCs, we can realize novel photonic components and devices, such as PhC waveguides, filters, demultiplexers, etc [1,2].An important condition for the application of PhC waveguides is the efficient light coupling between PhC waveguides and other conventional optical waveguides.For some PhC line-defect waveguides that comprise air-holes in a dielectric medium, good light coupling can be achieved by directly connecting the PhC waveguide with a dielectric slab waveguide [3].For PhC line-defect waveguides composed of dielectric rods, couplers with tapered structure, such as tapered PhC couplers [4,5] and tapered dielectric waveguides [6] can be used for efficient optical mode coupling.In this work, we consider the possibility of a PhC coupler that consists of varying-sized dielectric rods.Comparing with tapered PhC coupler that normally have 7∼10 layers of rods [4,5], the proposed PhC coupler is very compact, consisting of only 3 layers of dielectric rods along the light propagation direction.
The simulation and optimization of photonic components that contain both electrically large and small features is very challenging, because refined mesh grids around the small features may significantly increase the computational overhead [7].As shown in the design examples of section 3 and 4, PhC devices may contain dielectric rods that are much smaller than others, thus mesh refinements are needed to capture the local variations of the device geometry and electromagnetic field in conventional simulation and optimization process.This will lead to an increase of unknowns to be solved, and may cause problems of generating meshes of good quality around these small rods.The above difficulties, which are referred to as multi-scale problems [7], may become severe when the PhC device contains several dozens to hundreds of small rods.To avoid the multi-scale problems, we can apply the recently developed TO technique [7][8][9][10][11][12] to transform the original device model into an equivalent model, which consists of uniform and fixed-sized rods, with parameterized permittivity and permeability distribution.
Therefore, mesh refinement around small sized rods can be avoided, and PhC device can be simulated more efficiently.Photonic devices with varying geometry can be optimized with heuristic methods, such as the genetic algorithm [13] and the particle swarm algorithm [14].However, because of their high computational cost, heuristic methods are inefficient for devices with lots of design variables [15].In this work, we use a gradient method to optimize PhC devices, where the gradient of optimization object is calculated with the Adjoint Variable Method (AVM) [16][17][18].By using the AVM method, for each optimization iteration only one forward problem and one adjoint problem need to be solved.Therefore, the optimization is very efficient for devices with complex structures [15,16].

Transformed model of PhC device
Figure 1(a) shows the structure of a PhC device that consists of dielectric rods with varying size.To avoid the multi-scale problems, we use the TO technique [7][8][9][10][11] to transform the original device model into an equivalent model that consists of rods with a uniform and fixed size, as shown in Fig. 1(b).The transformed rods and surrounding medium have anisotropic permittivity and permeability distributions, parameterized by radii of the rods in real geometry.The parameterized permittivity and permeability distributions can be derived using the TO technique [7][8][9][10][11].Therefore, optimization of the original device is replaced by optimization of the parameterized permittivity and permeability distribution of the transformed device model.
A single dielectric rod and the transformed model is shown in Figs.1(c) and 1(d), respectively.The dielectric rod in real geometry has a radius of r c , and is stretched to be r b in the transformed model.The surrounding air of the rod in Fig. 1(c) is transformed into the surrounding medium layer of the stretched rod in Fig. 1(d), while the outer boundary of the surrounding medium is kept as same as r a after the transformation.Region of the real dielectric rod and surrounding air are denoted by Ω 1 and Ω 2 , respectively.In the transformed model, region of the stretched rod and the surrounding medium are denoted by Ω 1 and Ω 2 , respectively.Position vector in the original (transformed) geometry is denoted by r(x, y) (r( x, ỹ)).Therefore, the coordinate transformation from the original geometry of Fig. 1(c) to the transformed geometry of Fig. 1(d) can be defined as and With the above coordinate transformation, the permittivity and permeability distribution of the transformed model in Fig. 1(d) can be calculated as [7,11] εr where J is the Jacobian of the coordinate transformation In this work, we use commercial FEM software package COMSOL [19] to simulate PhC devices, with the transformed device model.To verify the transformed model, we first simulate light scattering from the original dielectric rod in Fig. 1(c), and then compare it with simulation using the transformed model in Fig. 1(d).In the simulation we set r a = 0.45µm, r b = 0.3µm, and r c = 0.125µm.Refractive index of the dielectric rod is n rod = 3.4.We assume the incident light propagates along x-axis, at the wavelength of λ = 1µm.The incident light polarizes along z-axis (TM mode), with the electric field intensity of |E| = 1.The total electric field obtained from simulations using the original and transformed model is shown in Figs.2(a) and 2(b), respectively.Comparing Fig. 2(b) with Fig. 2(a), we can see that the electric field distribution inside the transformed rod is stretched.However, in the region r > r a , the difference of the electric field between the two simulations is negligible (∼ 10 −4 ), as shown in Fig. 2(c).Thus for an outside observer, the transformed rod in Fig. 2(b) "looks" identical to the original rod in Fig. 2(a).
To evaluate the advantages of the TO based model, we compare the meshes of the original and transformed model of a single dielectric rod, in the region of r ≤ r b .r b is the radius of  we can see in order to accurately capture the small geometrical features, meshes around the rod in original model is much denser than the meshes of the transformed rod, especially for rod with smaller radius.Statistic data of meshes in region of r ≤ r b , including number of elements, mesh points, and degrees of freedom, are listed in Table 1.Data in Table 1 shows, for meshes of rod with radius of r c = 0.2 ∼ 0.001µm, the number of degrees of freedom is approximately 2 ∼ 7 times of that of the transformed model.From such a comparison, we can conclude that our TO based approach will be superior over conventional methods, especially for complex geometries with many variable and small features.The dispersion relation and guided mode of the slab waveguide can be obtained by solving the eigenvalue equation of the dielectric slab waveguide [20].

Optimization of PhC coupler
We assume the period of bulk PhC is a = 0.6µm.and the radius of the dielectric rod in bulk PhC is r 0 = 0.125a, with refractive index of n rod = 3.4.Band structure of the bulk PhC  is calculated using MIT Photonic Bands (MPB) [21], and Photonic Band Gap (PBG) is found in the range of normalized frequency 0.37335 ∼ 0.48741, which corresponds to wavelength 1.23µm ∼ 1.61µm.A PhC line defect waveguide is formed by removing one row of dielectric rods.Numerical simulations show that odd mode of the PhC line defect waveguide is supported in the wavelength range of 1.3µm ∼ 1.6µm.
To optimize PhC devices, an object function is defined as where p i is the power transmission coefficient of the device at the i-th sampling wavelength.
To evaluate the effectiveness of the proposed design method, we use three different sets of sampling wavelengths to define optimization objectives, namely ob j1 : λ i = 1.5, 1.525, 1.55, 1.575, 1.6(µm) ob j2 : λ i = 1.545, 1.55, 1.555(µm) (6) ob j3 : λ i = 1.55(µm)For sampling wavelengths set ob j1, which contains 5 different sampling wavelengths, the optimization procedure will check and optimize at all the 5 sampling wavelengths.Hence the optimized PhC device will have the broadest transmission band.For sampling wavelengths set ob j3, which contains only one sampling wavelength, the optimized PhC device will have the highest transmission coefficient at the sampling wavelength, but with the narrowest transmission band.
The optimization procedure is implemented by a free optimization package minConf [22], where the limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method is used to iteratively minimize the value of the object function J.For each optimization iteration, permittivity and permeability distribution of the transformed device model is calculated.Wave equations of the transformed model are solved to simulate light propagation in the device, and power transmission coefficient p i as well as the optimization object J are also calculated.Subsequently, values of r j for the next iteration are generated according to the first-order derivative term of the object ∂ J/∂ r j .The optimization iterations will exit when convergence of the object J is met.
The first-order derivative of the object function J subject to the design variable r j is calculated as where the design variable r j is the radius of the j-th rod of the PhC device.In the transformed device model r j represents a parameter of the permittivity and permeability distribution of the j-th transformed rod and its surrounding medium.The partial derivative term ∂ p i /∂ r j needs to be calculated for each design variable r j , and thus may results in intensive computational cost during the optimization procedure.In this work, the partial derivative term ∂ p i /∂ r j is calculated effectively by using the adjoint variable method [15][16][17][18], where merely one forward problem and one additional adjoint problem need to be solved to obtain the partial derivative term ∂ p i /∂ r j for all the design variables.To ensure the symmetry of the optimized PhC devices, the partial derivative term ∂ p i /∂ r j of dielectric rods at symmetric positions are averaged in the optimization process.
At the beginning of the optimization process, we set all the 18 rods in each coupler with radius of r 0 = 0.125a, which is the same as that of the bulk PhC.The transmission spectrum of the device with unoptimized PhC coupler is shown by the black dotted line in Fig. 5(a).Oscillations of the transmission spectrum of the unoptimized coupler indicates Fabry-Perot resonance caused by reflections at the unoptimized couplers, and the low power transmission coefficient (< 0.5) indicates an inefficient mode conversion by the unoptimized couplers.Figure 5(b) shows the profile of electric field E z of the unoptimized PhC coupler, calculated at the wavelength of 1.55µm.We can see clearly that light is radiated from the two unoptimized PhC couplers.
For the optimization with ob j1, ob j2, ob j3, the iteration step of optimization procedure is 71, 154, and 66, respectively.On a computer with Xeon E5645 CPU at 2.4GHz, the running time is 1016, 1302, 182 seconds, respectively.In this design example, at one sampling wavelength the computation time for the forward problem and the adjoint problem is approximately 2.3 and 0.5 seconds respectively.The optimized radii of the rods are shown in Fig. 4(b).Because of the symmetry of the couplers, only the radii of 9 rods of the coupler are listed in the table.The transmission spectra of the optimized coupler are shown in Fig. 5(a), where the spectra of couplers after optimization with ob j1, ob j2, and ob j3 are indicated by the dashed, solid, and dash-dotted line respectively.Figure 5(a) shows, the PhC coupler optimized with ob j1 has the broadest transmission band.For PhC coupler optimized with ob j3, the power transmission coefficient at the central wavelength of 1.55µm (specifically 0.9877) is slightly higher than the couplers optimized with ob j1 and ob j2 (0.945 and 0.9876 respectively).As a verification, we choose the coupler that optimized with ob j3, and simulate light propagation in the device.The obtained electric field E z is depicted in Fig. 5(c), clearly showing the almost perfect light coupling.

Optimization of PhC coupler and PhC waveguide bend
A schematic diagram of two PhC couplers together with a PhC waveguide bend is shown in Fig. 6(a).A PhC waveguide bend, indicated by dotted line, is located at the center of the device.Similar to the previous design example, we set the initial radii of the rods of the device to be r 0 = 0.125a.The same optimization objectives ob j1, ob j2, and ob j3 are used.For the optimization with ob j1, ob j2, ob j3, the iteration steps are 66, 236, and 149, with running time of 1544, 3275, 670 seconds, respectively.The optimized rods sizes are shown in Fig. 6(b).An interesting phenomenon is, for device that optimized with ob j3, the radius of rod of No. 12 declines to the preassigned lower bound of 0.001a (0.6nm).We can just neglect such a "tiny" dielectric rod, which is verified by further simulations.The spectra of the device are shown in Fig. 7(a).Similar to the previous design example, transmission band of the device that optimized with ob j1 is the broadest one, and power transmission coefficient of device that optimized with ob j3 is higher than others at the wavelength of 1.55µm.The electric field distribution E z of device that optimized with ob j3 are shown in Fig. 7(b).The rod of No. 12 and its counterpart at symmetric position are removed.

Perturbation analysis of the optimized device
The performance of PhC devices can be degraded by disorders originated from fabrication process [23,24].In this section, we will use sensitivity analysis to estimate the disturbance of the transmission of the optimized PhC device caused by fabrication imperfection of the rods sizes.To compare the disturbance of PhC devices optimized with different objectives, we define the normalized sensitivity as Δr(∂ p i /∂ r j ).If we interpret the normalization coefficient Δr as fluctuation of rods radii that originate from fabrication disorders, then the normalized sensitivity represents the disturbance of the power transmission coefficient p i , caused by the fabrication disorders.For the PhC devices optimized in previous sections, the magnitude of the normalized sensitivity are shown in Fig. 8, which are calculated at the wavelength of λ = 1.55µm.Symbols "+", "•", and "×" represent sensitivity for devices that are optimized with ob j1, ob j2, and  ob j3, respectively.The sensitivity is normalized with Δr = 0.01a, corresponding to 6nm of the fabrication error of the rods radii.Figure 8(a) shows, for the device optimized in section 3, the sensitivity of the device optimized with ob j1 is nearly two orders of magnitude larger than devices optimized with ob j2 and ob j3, thus it is more likely to be disturbed by fabrication disorders.At the same time, the rods of No. 1, 4, and 7, which are rods that close to the PhC line defect, are more sensitive than other rods.So more attention should be paid to these rods during the fabrication process.For PhC device optimized in section 4, the maximum magnitude of the normalized sensitivity of device optimized with ob j1, ob j2, and ob j3 is on the order of 10 −2 , 10 −3 , and 10 −4 respectively.For device that optimized with ob j1, the rods of No. 1, 4, 7 and 13 are the most sensitive rods of the device, as shown in Fig. 8(b).

Conclusions
We use the transformation optics modeling technique and gradient-based optimization method to simulate and optimize PhC devices with varying-sized dielectric rods.By using the transformation optics approach, the original device model is transformed into an equivalent model that consists of uniform and fixed-sized rods, with parameterized permittivity and permeability distributions.Therefore, mesh refinement around small rods can be avoided, and PhC devices can be simulated more efficiently.In addition, the adjoint variable method is used to calculated the gradient of the optimization object, which is very efficient for optimizing devices with multiple design variables.The proposed method is expected to find important applications in designing and optimizing other photonic and optical devices with complex structures.

Fig. 1 .
Fig. 1.(a) A PhC device consisting of varying-sized dielectric rods.(b) Transformed model of (a), in which the transformed rods have a uniform and fixed size.Gray color represents the permittivity distributions.(c) A single dielectric rod, surrounded by air.(d) Transformed model of (c).

Fig. 2 .
Fig. 2. (a) Light scattering of a single dielectric rod.(b) Light scattering of the same dielectric rod, using the transformed model.(c) The difference of the field distribution between (a) and (b) in the region of r > r a , showing almost negligible difference (in the order of 10 −4 ).

Figure 4 (
Figure4(a) is a schematic diagram of the proposed PhC coupler.An input coupler, denoted as coupler1, couples light between the input dielectric slab waveguide and the PhC line defect waveguide.Similarly, an output coupler, denoted as coupler2, couples light between the PhC waveguide and the output slab waveguide.The input and output coupler are symmetric, converting the slab waveguide mode into the PhC line defect waveguide mode or vice versa.Both the input and output slab waveguide are single mode waveguides, with core thickness of d = 1µm, and refractive index of n co = 1.5.The slab waveguide core is surrounded by air.The dispersion relation and guided mode of the slab waveguide can be obtained by solving the eigenvalue equation of the dielectric slab waveguide[20].We assume the period of bulk PhC is a = 0.6µm.and the radius of the dielectric rod in bulk PhC is r 0 = 0.125a, with refractive index of n rod = 3.4.Band structure of the bulk PhC

Fig. 4 .
Fig. 4. (a) Schematic diagram for PhC coupler.PhC couplers, indicated by dashed lines, couple light between the PhC line defect waveguide and the dielectric slab waveguides.(b) Rods radii of the optimised coupler, which are optimized with sampling wavelengths ob j1-ob j3, as defined in Eq. (6).

Fig. 5 .
Fig. 5. (a) Transmission spectra of the PhC coupler for a straight PhC waveguide.Spectrum of the unoptimized coupler is indicated by dotted line (black), spectrum of the coupler that optimized with ob j1, ob j2, and ob j3 is indicated by dashed line (green), solid line (blue), and dash-dotted line (red) respectively.(b) Electric field E z of the unoptimized PhC coupler, calculated at the wavelength of 1.55µm.(c) Electric field E z of the coupler optimized with ob j3.

Fig. 6 .Fig. 7 .Fig. 8 .
Fig. 6.(a) Schematic diagram for PhC coupler and PhC waveguide bend.PhC couplers are indicated by dashed lines, and PhC waveguide bend is indicated by dotted line.(b) Rods radii of the optimised PhC coupler and PhC waveguide bend, which are optimized with sampling wavelengths ob j1-ob j3, as defined in Eq. (6).

Table 1 .
Statistic data of meshes of the original and transformed model, in the region of r ≤ r b .Unit for radius: µm.