Inverse designed photonic crystal de-multiplex waveguide coupler

A two-dimensional photonic crystal diplexer integrated with a waveguide coupler is proposed. The design is computer generated through an inverse design process, limited within an area measuring 5μm× 5μm. The best working device was designed for the optical communication wavelengths, 1.50μm and 1.55μm, i.e. a channel spacing of 50 nm. The device exhibits crosstalks suppressed below 40dB and coupling efficiencies close to 80%, for both channels. © 2005 Optical Society of America OCIS codes: (130.2790) Guided waves; (130.3120) Integrated optics devices; (060.1810) Couplers, switches, and multiplexers; (220.4830) Optical systems design; (999.9999) Photonic Crystals References and links 1. E. Yablonovitch, ”Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. 58, 2059-2062 (1987). 2. S. John, ”Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 24862489 (1987). 3. J. D. Joannopoulos, R. D. Meade and J. N. Winn, Photonic Crystals, (Princeton Press, Princeton, New Jersey, 1995). 4. L. Sanchis, A. Håkansson, D. Lopez-Zanón, J. Bravo-Abad and J. Sánchez-Dehesa, ”Integrated optical devices design by genetic algorithm,” Appl. Phys. Lett. 84, 4460–4462 (2004). 5. L. Shen, Z. Ye, and S. He, ”Design of two-dimensional photonic crystals with large absolute band gaps using a genetic algorithm,” Phys. Rev. B 68, no. 035109 (2003). 6. Stefan Preble, Michal Lipson and Hod Lipson, ”Two-dimensional photonic crystals designed by evolutionary algorithms,” Appl. Phys. Lett. 86, 061111-061113 (2005). 7. Burger M., S. J. Osher, and E. Yablonovitch, ”Inverse problem techniques for the design of photonic crystals,” IEICE Trans. Electron. E87C, 258–265 (2004). 8. J. Geremia, J. Williams, and H. Mabuchi, ”Inverse-problem approach to designing photonic crystals for cavity QED experiments,” Phys. Rev. E 66, no. 066606 (2002). 9. I. Gheorma, S. Haas, and A. Levi, ”Aperiodic nano-photonic design,” J. of Appl. Phys. 95, 1420-1426 (2004). 10. P. Borel, A. Harpoth, L. Frandsen, M. Kristensen, P. Shi, J. S. Jensen and O. Sigmund, ”Topology optimization and fabrication of photonic crystal structures,” Opt. Express 12, 1996-2001 (2004). 11. Jasmin Smajic, Christian Hafner and Daniel Erni, ”Optimization of photonic crystal structures,” J. Opt. Soc. Am. A 21, 2223–2232 (2004). 12. T.D. Happ, M. Kamp and A. Forchel, ”PhC tapers for ultracompact mode conversion,” Opt. Lett. 26, 1102–1104 (2001). 13. P. Bienstman , S. Assefa, S.G. Johnson, J. D. Joannopoulos, G. S. Petrich and L. A. Kolodziejski , ”Taper structures for coupling into PhC slab waveguides,” J. Opt. Soc. Am. B 20, 1817–1821 (2003). 14. P. Sanchis, J. Mart, A. Garca, A. Martnez and J. Blasco, ”High efficiency coupling technique for planar photonic crystal waveguides,” Electron. Lett. 38, 961–962 (2002). 15. S. Noda, A. Chutinan and M. Imada, ”Trapping and emission of photons by a single defect in a photonic bandgap structure,” Nature (London) 407, 608–610 (2000). 16. H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Tamamura, T. Sato and Kawakami S, ”Photonic crystals for micro lightwave circuits using wavelength-dependent angular beam steering,” Appl. Phys. Lett. 74, 1370–1372 (1999). (C) 2005 OSA 11 July 2005 / Vol. 13, No. 14 / OPTICS EXPRESS 5440 #7504 $15.00 USD Received 20 May 2005; revised 28 June 2005; accepted 28 June 2005 17. A. Håkansson and José Sánchez-Dehesa, ”Inverse design of photonic crystal devices,” IEEE J. Sel. Area Comm. (2005) (to be published). 18. A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering, (Prentice Hall, New Jersey, 1991). 19. Davy Pissoort, Bart Denecher, Peter Bienstman, Fank Olyslager and Danil De Zutter, ”Comparative study of three methods for the simulation of two-dimensional photonic crystals,” J. Opt. Soc. Am. A 21, 2186–2195 (2004). 20. Davy Pissoort and Frank Olyslager, ”Termination of periodic waveguides by PMLs in time-harmonic integral equation-like techniques,” IEEE Antennas and Wireless Propagation Letters 2, 281–284 (2003). 21. D.E. Goldberg, Genetic Algorithms in Search, Optimization and Learning, (Addison Wesley, Reading , MA, 1989). 22. J.H. Holland, Adaptation in natural and Artificial Systems, (The University of Michigan Press, Ann Arbor, 1975). 23. B.R. Moon, Y.S. Lee and C.K. Kim, ”GEORG: VLSI circuit partitioner with a new genetic algorithm framework,” Journal of Intelligent Manufacturing 9, 401-412 (1998). 24. Andreas Håkansson, José Sánchez-Dehesa, ”Comment on ’Optimization of photonic crystal structures’, J. Opt. Soc. Am. A 21, 2223-2232 (2004),” arXiv.org, cond-mat/0504581 (2005) 25. E. Cantú-Paz, D. E. Goldberg, ”Are multiple runs of genetic algorithms better than one?,” Genetic and Evolutionary Computation GECCO 2003, PT I, Proceedings lecture notes in Computer Science 2723, 801–812 (2003). 26. S. Boscolo and M. Midrio, ”Three-dimensional multiple-scattering technique for the analysis of photonic-crystal slabs,” J. Lightwave Tech. 22, 2778–2786 (2004).


Introduction
Photonic crystals (PhCs) [1,2,3] are new materials with unprecedented properties for microscaled optical device design.By introducing defects in these crystalline materials light can be manipulated and guided through sub-wavelength cavities with very low losses.These new components are normally designed by physical intuition, based on the insight of the nature of PhCs.This direct design method, from property to functionality, guided by intuition, is limited in many ways.To reach beyond our knowledge of design, to find unexplored solutions, the guidance has to be done, not by intuition, but by methodically optimization.Instead of looking for functionalities from a set of properties, we search for the device for a fixed functionality.The inverse design approach is based on defining the wanted functionality a priori, then, a functional device is optimized within set constrains.In the field of PhCs applications, ID has been used to successfully design, e.g spot-size converters [4], photonic band-gap materials [5,6,7], cavities for QED-experiments [8], aperiodic scatterers for tailored transmission or reflection [9], low loss PhC wave-guide Z-bend [10] and PhC power splitters [11].
In this work we apply our ID method, introduced in [4], to design an ultra compact PhC 1 × 2 de-multiplex coupler for separating and coupling two wavelength from a single dielectric WG to two separate PhC-WGs.Each of these two functionalities have earlier been solved separately using the properties of PhCs.The more promising direct designed proposals for coupling a dielectric WG to a PhC-WG are the use of tapers [12,13,14].Moreover, the manipulation of the light propagating inside the PhC structure has shown fruitfulness regarding micro scaled wavelength-division-multiplexing (WDM).By integrating isolated defects with the PhC-WG, Noda et al. demonstrated the feasibility of compact add/drop filters using the trapping and emission of photons [15].A second approach for WDM by using light-beam steering was proposed by Kosaka et al. [16].
We will in this work show that it is possible to design a device that solves both of these two functionalities without any extensive knowledge about PhC properties.The features of this 'two devices in one' design can be compared with the best proposed direct designed PhC-WG couplers as well as the WDM devices.

Inverse design method
A fictional state-of-the-art ID-process would consist in, one input that defines the available resources, such as, fabrication and financial limitations, and a second input defining the askedfor functionality.Using these two inputs the process simply finds the best possible design with the set functionality, through rigor optimization and simulations.Unfortunately this includes too complex calculations, and the problem has to be limited and scaled down to be practical and possible to simulate.
In the following two subsection our choice of simulation and optimization method will be described.A more detailed survey of this ID approach can be found in [17].

Simulation
Multiple scattering theory (MST) [18] is one of the common methods for wave scattering simulations.This semi-analytic method has various advantages [19] over more popular wave simulation tools, such as the finite-difference time-domain method (FDTD).The MST consists in expanding the propagating fields, induced by each scatterer, in a base of natural Bessel-functions.Using a Bessel functions expansion gives the advantage of a single parameter control of the accuracy, as well as a fast convergence for circular scattering objects.This fast convergence leads to very short simulation times for clusters of a limited numbers of scattering objects.Though the MST is considered fast, the computational complexity scales as the cube of the number of cylinders, which implies that this parameter always has to be kept under consideration.Furthermore, for a finite number of scatterers the method does not rely on an absorbing boundary conditions.Consequentially, the theory is suited, though not limited, for clusters with a fixed number of scatterers.Although, infinite, or semi-infinite, structures can be simulated by introducing materials with absorption.Pissoort et al. [20] showed how a PhC-based perfectly matched layer (PML) can be implemented by introducing an imaginary part in scatterer's position-coordinates.
Another direct solvers used for ID processing PhCs is the frequency domain finite element model.The disadvantages for using numerical methods are the difficulties of error control introduced by the phase, as well as the additional memory and CPU costs.For the structures treated here the MST approach is superior in both accuracy and speed.On the other hand numerical methods do not depend on the shape of each scatterer and are well suited for optimizing arbitrary shapes, e.g.see [10].

Optimization
Basically, any optimization method could be used to solve the ID optimization problem, though the choice will set demands on the direct solver as well as restrictions on the design variable.The optimization algorithm used in this ID-process is a Genetic Algorithm (GA) [21,22].GAs are very popular in engineering optimization problem, in part for their superior performance and in part for their straight forward implementation.Even though GAs have proved to be very competent in finding global maximums for complex optimization problems, they normally require many evaluations before convergence is reached.This will be a problem if the time needed to evaluate the optimization function (fitness) of a design is not short enough.In the case of MST, the simulation time is proportional to N 3 , where N is the number of cylinders.Hence, a suitable speed can be achieved if the number of cylinders is kept sufficiently small.To be able to run the PhC structure optimization on a standard Pentium-IV processor the maximum number of cylinders has here been set to approximately 200-250.A second restriction on the design variable, set by the optimization process, is the discrete search space.GAs use binary implementation of the optimization parameters, which practically excludes float parameters.
Here the binary parameters represent the absence or presence of different cylinders.
The functionality of a GAs is easily explained by using its mimic of the natural selection of the evolution.It works simultaneously with a set of solutions ('population'), where each solution is coded as a 'genome' on a virtual chromosome.This genome is normally fed to a computer as a binary string of bits.In the setup for the PhC structure problem, each bit ('gene') on the chromosome correspond to a fixed lattice site and the binary value of the gene ('allele') represents the presents or the absence of a cylinder at that fixed position.The optimization normally takes off from a random initiation, i.e. there is no need for an initial starting point.The 'evolution' is then guided from generation to generation by three operators; Selection, Crossover and Mutation.The Selection simply selects better solution over worse, for 'mating'.The mating act is directed by the Crossover that mixes the alleles from two parent genomes, creating an offspring.Before placing the offspring in the new generation's population, the Mutation changes randomly a small percentage of the alleles on the chromosome.The advance in generations is complete when the new population is equal in size as the initial one.b) The repetition of the lattice is marked out with squares, representing the lattice sites.c) A two dimensional representation of the virtual genome.Each lattice site is filled with a binary digit, where "0" and "1" represents the presence or absence of a cylinder, respectively.d) The color scale represents the geographical linkage with respect to the central allele.Blue (red) corresponds to high (low) geographical linkage.
Intuitively, lattice sites placed close to one another should have higher correlation, i.e. higher geographical linkage [23], than lattice sites placed farther apart.Hence, if coding the chromosome by using a direct 2D representation as illustrated in Fig. 1, a two-dimensional crossover could be used efficiently.Here the 2D-geographic crossover [23] is used.The crossover simply works by slicing up the 2D-chromosome into clusters of neighboring cylinders.These clusters mark the boundaries of the interchange between the two parents.

Design variables
The design variables in the case of optimizing a de-multiplex waveguide coupler (DEMUX-WGC) are set by the coupling efficiency and the crosstalk attenuation, for respective channel and wavelength.To estimate these parameters, in a simple and accurate way, the amplitudes of the electric field in the center of the PhC-WGs are calculated.Since the parameters for the PhC-WG are chosen to obtain a single mode WG (TM-modes), the maximal amplitude of the mode-profile can be scaled proportionally with the total intensity flux of the coupled mode.This will result in four electric fields modulus (out-of-plane oscillations), E 11 , E 12 , E 21 and E 22 , where the first index correspond to λ 1 or λ 2 , and the second to channel-1 or channel-2.Using these four variables the coupling efficiency, CE, and the cross talk attenuation, XT , are calculated as, The n 1 and n 2 indexes each corresponds to one of the two channels.The reason for not fixing a specific wavelength to a fixed waveguide is to help the GA at its initial generations to faster find a promising solution.If after the ID process the wavelengths are coupled to the wrong channels, one can simple flip the structure with respect to the x-axis to solve this problem.However, if the first generations of the optimization is dominated by two promising solutions with inverse n 1 and n 2 indexes, the process could momentarily get stuck, before one solution gets dominant.In this work, nevertheless, this short delay has not been a significant problem.Though, if the consequences are increased, a simple inversion function could easily be implemented to operate on all "inverted" solutions by flipping the structure with respect to the x-axis, before calculating CE and XT .In this way all solutions will be directly comparable and n 1 and n 2 can each be fixed to an individual channel.
The fitness function has to maximize the coupling efficiency as well as the crosstalk for both wavelengths.A naive linear implementation of this fitness would probably lead to undesired results, where one of the two parameters gets dominant.In this case the crosstalk has a tendency to dominate over the coupling efficiency.To hamper this dominant effect it is necessary to set up a cut-off limit at a good-enough value, for which higher values of the cross talk do not increase significantly the total fitness.This can be implemented through nonlinear fitness scaling.In this work a fitness scaling of the second order is used for the cross talk parameters and has been implemented as, where α controlls the value of the cut-off value.Please notice that XT in Eqs. ( 1) and ( 2) correspond to the cross talk attenuation and consequentially the parameter has to be maximized and not minimized.The fitness parameter used in the optimization is now defined as, Figure 2 prints out the fitness dependency for one of the two addends in Eq. ( 3).The nonlinear curvature of the surface, with respect to the cross talk variable, can clearly be observed.

Optimization setup
Since we are trying to solve the inverse problem of coupling an external beam into a semiinfinite PhC-WG, a perfectly matching layers (PML) has to be implemented to run accurate simulations using the MST.Due to the limit of 250 cylinders, it is very important to minimize the total number of cylinders in the PML structure, which has to be implemented with as few Fig. 2. The fitness as a function of the coupling efficiency and the crosstalk.The nonlinear scaleing of the crosstalk parameter is implemented as Eq. ( 2) with α = 5. cylinders as possible.Pissoordt et.al [20] showed that if the complex plane of the PML is entered in a smooth way the reflection at the interface is strongly reduced and the wanted "semiinfinite" effect can be achieved with a much faster convergence.They proposed a circular PML where the imaginary part increases as "1 − cos".We minimized the size of this circular PML, adjusted for our PhC-WG, to a length of four lattice constants.The complex lattice points for this minimized PML, used in the positive x-direction, starting at x 0 = 0.0 and for a hexagonal lattice PhC-WG, with a missing line defect, was implemented as, were a is the lattice constant and n = 1, 2,...,8 is the file index.
The two PMLs, one for each PhC-WG, result in a total of 64 cylinders.In the setup the waveguides are positioned with a 60 degrees separation (±30 degrees with respect to the incident wave) to avoid coupling between the WG-modes.An extra 47 cylinders have been added to the WGs in order to assure that the electric field inside the PhC-WGs actually correspond to a proper propagating WG-mode.These 111 cylinders are set as the fixed structure in our simulations and remain intact throughout the optimization process.
The spacial distribution of the lattice sites, which to include in the optimization process, can be chosen in many ways and can be quite difficult to assert intuitively.In particular, the lattice sites need to be able to resemble two functionalities through the optimization.First, to focus; the incident beam has to be stretched in order to be coupled into the sub-wavelength PhC-WGs.Secondly to diplex; the light has to be coupled through cavities with high enough Q-value to be able to separate the two working wavelengths.The initial approach was done by (C) 2005 OSA 11 July 2005 / Vol. 13, No. 14 / OPTICS EXPRESS 5445 using two separated clusters, each cluster chosen to solve each functionality.A fixed numbers of files were set to the left hand side of the PhC-WGs, wide enough to be able to manipulate the total incident energy, with the intention to overcome the focusing problem, see Cluster-1 in Fig. 3. To minimize the reflectance, the lattice constant for this cluster has been modified so the correspondent infinite crystal has a propagating mode for the incident direction [4].A second cluster, Cluster-2, was placed at the opening of the two waveguides making it possible to form coupling-by-modes.Here the internal crystal structure is simply an expansion of WG-Fig.3. The first crystal structure optimization set-up.The blue dots marks the cylinders implemented with absorption for semi-infinite WG simulation.The green dots correspond to the PhC-WG and are fixed throughout the optimization process.Cluster-1 and Cluster-2, identified by the black and red dots/circles are the lattice sites represented by the virtual genome in the GA.The optimized structure is identified by the dots (alleles=1), the circles correspond to absent lattice sites (alleles=0).The scale used is normalized to the lattice parameter a.
crystal structure.A similar setup has successfully been applied in designing a PhC-WG coupler [4].The full lattice position setup is displayed in Fig. 3 with the optimized crystal structure.The resulting device from the ID-procedure gives an indication that this approach is not optimal.Though the reasoning for Cluster-2 seems to be based on accurate assumptions, Cluster-1 did not.The best functional device for this configuration includes an aperture at x = 0.0 in Cluster-1, as demonstrated in the figure.The fact that the most important lattice sites, at the high intensity area for the incident beam, are simply left vacant concludes that these positions do not reflect much importance in the designing process and, consequentially, the positions in Cluster-1 should be reconsidered.The modification made for the second approach are based on these aforementioned observations.The number of lattice sites placed in the mouth of the PhC-WG has been increased in number so the total incident beam can now be manipulated by this cluster.Since

Results
The calculations were done considering a hexagonal PhCs structures based on silicon rods (n Si = 3.4) placed in a background of silica (n SiO 2 = 1.46).The radius of the rods was set to 0.2a, where a is the lattice constant.In this crystal, the band gap for TM-modes opens for the frequency range of 0.266 − 0.359(a/λ ), where λ is the wavelength in free space.Moreover, the single line-defect PhC-WG has a guided single-mode for the frequency range 0.274 − 0.347(a/λ ).The width of the incident beam that is to be coupled into the PhC-WG was set to 15 lattice-constants.If a 0 = 0.465μm, this equals ≈ 6μm corresponding to the width of a dielectric WG.The two working-wavelengths were set, with a 50nm channel spacing, to λ 1 = 1.50μm and λ 2 = 1.55μm, corresponding to the guiding frequencies a 0 /λ 1 = 0.31 and a 0 /λ 2 = 0.30, respectively.Figure 5 shows two graphs that illustrate the best ID result for the lattice sites distribution setup illustrated in Fig. 4 and for the optimization function in Eq. (3) using the scaling constant α = 5.The GA optimization was carried out using a population of 200 genomes, the Tournament-3 Selection [21], the 2d-geographical Crossover [23] and the Mutation set to 0.66%.Regarding the population size used in the optimization, multiple runs using a small population can be perferable in less complex PhC structure optimization [11], however when dealing with larger genomes a larger population is in general superiour [24,25].The GA-run converged after 50'000 fitness evaluations (250 generations), which equals an effective CPUtime of approximately 120h on a 3.0GHz Pentium-IV processor.Fig. 5.The amplitude of the electric field for (left) a 0 /λ 1 = 0.31 and (right) a 0 /λ 1 = 0.30, where red (blue) corresponds to a large (small) amplitude.The scale used is normalized to the lattice parameter a.
To measure the transmission characteristics of the optimized PhC DEMUX-WGC, the power flux through a 2 √ 3a long line segment, placed at the center of the PhC-WG, was calculated and compared to the total incident power.The resulting coupling efficiencies and cross-talks for both wavelengths are printed out in Table 1.
Figure 6 shows the coupling into the two wave guide ports and reflection for the guided frequency spectra of the PhC-WG.The values are normalized to the total incident power through the input port.The two working wavelengths, marked by 0.30a/λ and 0.31a/λ in the figure are clearly marked out with a maximum and minimum, as expected for an optimized structure.One can also observe that the device experience almost total reflection for the frequencies outside the working wavelengths.

Conclusion
We have presented an optimized PhC device that, first, couples an external dielectric WG to two PhC-WGs, and secondly, diplex two wavelengths by coupling each wavelength to a different channel.This ultra compact device, only 5 μm wide, shows great features for optical wavelength communication.The predicted crosstalk was pressed below 40dB for both channels with a coupling efficiency close to 80%.Although these simulations was done in a 2D sub-space, these results confirm a necessary and important step toward full 3D-ID.It is clear that fabricated 2D-PhC slab structures show, to some extent, similarities with 2D simulations.This fact Fig. 6.Normalized intensity spectra of the DEMUX-WGC device in Fig. 4. The blue line corresponds to port 1 (the top wave guide), the red line to port 2 (the bottom wave guide) and the black line show the loss due to reflection.The power is normalized to the total incident power through the input port.
states that a 3D-ID process should always be based on a initial 2D-ID result.The MST can be extended to full 3D simulation for PhC-slabs [26] by amplify the field expansions to include the field profiles along the z-axis.This will considerable increase the demand on computational power, but if a 2D-ID is carried out a priori, a 3D-ID process is realistic within a reasonable timeframe.

Fig. 1 .
Fig. 1.GA representation of a PhC structure.a) The black dots correspond to the PhC structure.b)The repetition of the lattice is marked out with squares, representing the lattice sites.c) A two dimensional representation of the virtual genome.Each lattice site is filled with a binary digit, where "0" and "1" represents the presence or absence of a cylinder, respectively.d) The color scale represents the geographical linkage with respect to the central allele.Blue (red) corresponds to high (low) geographical linkage.

Fig. 4 .
Fig.4.The second crystal structure optimization set-up.The blue dots marks the cylinders implemented with absorption for semi-infinite WG simulation.The green dots correspond to the PhC-WG and are fixed throughout the optimization process.Cluster-1, identified by the black dots/circles are the lattice sites represented by the virtual genome in the GA.The optimized structure is identified by the dots (alleles=1), the circles correspond to absent lattice sites (alleles=0).The scale used is normalized to the lattice parameter a.
this cluster (C) 2005 OSA 11 July 2005 / Vol. 13, No. 14 / OPTICS EXPRESS 5446 now is wide enough, it can also function in a taper-like way; the focusing and the coupling problems can both be solved by Cluster-1, see Fig 4. The new optimizations resulted in a increase of 138% of the optimized fitness value.