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
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 . For PhC line-defect waveguides composed of dielectric rods, couplers with tapered structure, such as tapered PhC couplers [4, 5] and tapered dielectric waveguides  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 . 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 , 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–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  and the particle swarm algorithm . However, because of their high computational cost, heuristic methods are inefficient for devices with lots of design variables . 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–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].
2. 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–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–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 rc, and is stretched to be rb 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 ra 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 asFig. 1(d) can be calculated as [7, 11]
In this work, we use commercial FEM software package COMSOL  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 ra = 0.45μm, rb = 0.3μm, and rc = 0.125μm. Refractive index of the dielectric rod is nrod = 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 > ra, 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 ≤ rb. rb is the radius of the transformed rod, which is set as rb = 0.3μm as above. Meshes in the region of r ≥ rb can be generated identically for the original and transformed model, so they are ignored in the comparison. In the region of r ≤ rb, triangular mesh grids are generated, and quadratic elements are used. Meshes of the transformed rod, which just occupy the region of r ≤ rb, are shown in Fig. 3(a). In the region of r ≤ rb, meshes of the original model are shown in Figs. 3(b)–3(e), with rod radius of rc = 0.2, 0.1, 0.05, 0.01μm, respectively. Comparing Fig. 3(a) with Figs. 3(b)–3(e), we can see that 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 ≤ rb, 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 rc = 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.
3. Optimization of PhC coupler
Figure 4(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 nco = 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 .
We assume the period of bulk PhC is a = 0.6μm. and the radius of the dielectric rod in bulk PhC is r0 = 0.125a, with refractive index of nrod = 3.4. Band structure of the bulk PhC is calculated using MIT Photonic Bands (MPB) , 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
The optimization procedure is implemented by a free optimization package minConf , 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 pi as well as the optimization object J are also calculated. Subsequently, values of rj for the next iteration are generated according to the first-order derivative term of the object ∂J/∂rj. 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 rj is calculated as15–18], where merely one forward problem and one additional adjoint problem need to be solved to obtain the partial derivative term ∂pi/∂rj for all the design variables. To ensure the symmetry of the optimized PhC devices, the partial derivative term ∂pi/∂rj 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 r0 = 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 Ez 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 obj1, obj2, obj3, 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 obj1, obj2, and obj3 are indicated by the dashed, solid, and dash-dotted line respectively. Figure 5(a) shows, the PhC coupler optimized with obj1 has the broadest transmission band. For PhC coupler optimized with obj3, the power transmission coefficient at the central wavelength of 1.55μm (specifically 0.9877) is slightly higher than the couplers optimized with obj1 and obj2 (0.945 and 0.9876 respectively). As a verification, we choose the coupler that optimized with obj3, and simulate light propagation in the device. The obtained electric field Ez is depicted in Fig. 5(c), clearly showing the almost perfect light coupling.
4. 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 r0 = 0.125a. The same optimization objectives obj1, obj2, and obj3 are used. For the optimization with obj1, obj2, obj3, 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 obj3, 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 obj1 is the broadest one, and power transmission coefficient of device that optimized with obj3 is higher than others at the wavelength of 1.55μm. The electric field distribution Ez of device that optimized with obj3 are shown in Fig. 7(b). The rod of No. 12 and its counterpart at symmetric position are removed.
5. 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(∂pi/∂rj). 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 pi, 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 obj1, obj2, and obj3, 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 obj1 is nearly two orders of magnitude larger than devices optimized with obj2 and obj3, 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 obj1, obj2, and obj3 is on the order of 10−2, 10−3, and 10−4 respectively. For device that optimized with obj1, the rods of No. 1, 4, 7 and 13 are the most sensitive rods of the device, as shown in Fig. 8(b).
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.
This work is supported by the National Natural Science Foundation of China (Grand No. 51275504), and the Startup Fund for Postdoctoral Scientific Research Program of Jilin Province.
References and links
1. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2 (Princeton University, 2008).
2. I. A. Sukhoivanov and I. V. Guryev, Photonic Crystals: Physics and Practical Modeling (Springer-Verlag, Berlin, 2009). [CrossRef]
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). [CrossRef]
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). [CrossRef]
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). [CrossRef]
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). [CrossRef]
7. O. Ozgun and M. Kuzuoglu, “Software metamaterials: Transformation media based multiscale techniques for computational electromagnetics”, J. Comput. Phys. 236, 203–219 (2013). [CrossRef]
8. A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equation,” J. Mod. Opt. 43(4), 773–793 (1996). [CrossRef]
11. L. H. Gabrielli, D. Liu, S. G. Johnson, and M. Lipson, “On-chip transformation optics for multimode waveguide bends,” Nat. Commun. 3:1217(2012). [PubMed]
12. Q. Wu, J. P. Turpin, and D. H. Werner, “Integrated photonic systems based on transformation optics enabled gradient index devices,” Light: Science & Applications 1, e38; (2012). [CrossRef]
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). [CrossRef] [PubMed]
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). [CrossRef] [PubMed]
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). [CrossRef]
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). [CrossRef]
19. COMSOL Multiphysics, http://www.comsol.com/
20. D. Marcuse, Theory of Dielectric Optical Waveguides(Academic Press, 1974).
22. M. Schmidt, http://www.di.ens.fr/~Emschmidt/Software/minConf.html
24. V. Savona, “Electromagnetic modes of a disordered photonic crystal,” Phys. Rev. B 83, 085301 (2011). [CrossRef]