## Abstract

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

## 1. 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–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–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 *r _{c}*, and is stretched to be

*r*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

_{b}*r*after the transformation. Region of the real dielectric rod and surrounding air are denoted by Ω

_{a}_{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

**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*= 0.3μm, and

_{b}*r*= 0.125μm. Refractive index of the dielectric rod is

_{c}*n*= 3.4. We assume the incident light propagates along

_{rod}*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*, the difference of the electric field between the two simulations is negligible (∼ 10

_{a}^{−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*is the radius of the transformed rod, which is set as

_{b}*r*= 0.3μm as above. Meshes in the region of

_{b}*r*≥

*r*can be generated identically for the original and transformed model, so they are ignored in the comparison. In the region of

_{b}*r*≤

*r*, triangular mesh grids are generated, and quadratic elements are used. Meshes of the transformed rod, which just occupy the region of

_{b}*r*≤

*r*, are shown in Fig. 3(a). In the region of

_{b}*r ≤ r*, meshes of the original model are shown in Figs. 3(b)–3(e), with rod radius of

_{b}*r*= 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

_{c}*r*≤

*r*, 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

_{b}*r*= 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.

_{c}## 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 *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.125*a*, 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*is the power transmission coefficient of the device at the

_{i}*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

*obj*1, 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

*obj*3, 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*for the next iteration are generated according to the first-order derivative term of the object

_{j}*∂J/∂r*. The optimization iterations will exit when convergence of the object

_{j}*J*is met.

The first-order derivative of the object function *J* subject to the
design variable *r _{j}* is calculated as

*r*is the radius of the

_{j}*j*-th rod of the PhC device. In the transformed device model

*r*represents a parameter of the permittivity and permeability distribution of the

_{j}*j*-th transformed rod and its surrounding medium. The partial derivative term

*∂p*needs to be calculated for each design variable

_{i}/∂r_{j}*r*, and thus may results in intensive computational cost during the optimization procedure. In this work, the partial derivative term

_{j}*∂p*is calculated effectively by using the adjoint variable method [15–18], where merely one forward problem and one additional adjoint problem need to be solved to obtain the partial derivative term

_{i}/∂r_{j}*∂p*for all the design variables. To ensure the symmetry of the optimized PhC devices, the partial derivative term

_{i}/∂r_{j}*∂p*of dielectric rods at symmetric positions are averaged in the optimization process.

_{i}/∂r_{j}At the beginning of the optimization process, we set all the 18 rods in each coupler
with radius of *r*_{0} = 0.125*a*,
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 *obj*1, *obj*2,
*obj*3, 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 *obj*1, *obj*2, and
*obj*3 are indicated by the dashed, solid, and dash-dotted line
respectively. Figure 5(a) shows, the PhC
coupler optimized with *obj*1 has the broadest transmission band. For
PhC coupler optimized with *obj*3, the power transmission coefficient
at the central wavelength of 1.55μm (specifically 0.9877) is slightly higher
than the couplers optimized with *obj*1 and *obj*2
(0.945 and 0.9876 respectively). As a verification, we choose the coupler that
optimized with *obj*3, 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.

## 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
*r*_{0} = 0.125*a*. The same
optimization objectives *obj*1, *obj*2, and
*obj*3 are used. For the optimization with *obj*1,
*obj*2, *obj*3, 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 *obj*3, the
radius of rod of No. 12 declines to the preassigned lower bound of
0.001*a* (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 *obj*1 is the broadest one, and power transmission
coefficient of device that optimized with *obj*3 is higher than
others at the wavelength of 1.55μm. The electric field distribution
**E _{z}** of device that optimized with

*obj*3 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*(*∂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*, 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

_{i}*λ*= 1.55μm. Symbols “+”, “○”, and “×” represent sensitivity for devices that are optimized with

*obj*1,

*obj*2, and

*obj*3, respectively. The sensitivity is normalized with Δ

*r*= 0.01

*a*, 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

*obj*1 is nearly two orders of magnitude larger than devices optimized with

*obj*2 and

*obj*3, 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

*obj*1,

*obj*2, and

*obj*3 is on the order of 10

^{−2}, 10

^{−3}, and 10

^{−4}respectively. For device that optimized with

*obj*1, the rods of No. 1, 4, 7 and 13 are the most sensitive rods of the device, as shown in Fig. 8(b).

## 6. 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.

## Acknowledgments

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]

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

**10. **Y. M. Liu and X. Zhang, “Recent advances in transformation
optics,” Nanoscale **4**(17),
5277–5292 (2012). [CrossRef] [PubMed]

**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]

**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). [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]

**17. **G. Veronis, R. W. Dutton, and S. Fan, “Method for sensitivity analysis of
photonic crystal devices,” Opt.
Lett. **29**(19),
2288–2290 (2004). [CrossRef] [PubMed]

**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).

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**24. **V. Savona, “Electromagnetic modes of a disordered
photonic crystal,” Phys. Rev. B **83**, 085301 (2011). [CrossRef]