## Abstract

In this work, we use the finite-difference time-domain method in conjunction with a genetic algorithm (GA) to design a photonic crystal waveguide with gratinglike surface added for highly-efficient directional emission. By placing a detector at different locations in the output field, we have obtained both on-axis and off-axis highly-efficient directional emission designs with the help of GA. The interference of light emitted from the outlet of waveguide and surfaces modes in gratinglike surface is believed to account for the directional beaming effect in our designs.

©2009 Optical Society of America

## 1. Introduction

The diffraction limit is a basic principle in classical optics. As we know, because of the existence of diffraction limit, the emission of light out of a structure with a subwavelength feature size will spread into a wide angle range, and this is an undesirable property in near-field optics and integrated optics. Recently, both theoretical [1] and experimental [2] efforts have demonstrated that on-axis directional emission of light can be realized via a photonic crystal waveguide (PCW), if the surface of the waveguide at the output side, called endface or termination, is modified by a periodic corrugation. This directional emission is mainly attributed to the excitation of surface modes, which can couple to the continuum of air modes. However, due to reflection and impedance mismatch at the termination of PCWs, the reported directional emissions are inefficient. Therefore, many succedent researchers [3–11] focused there efforts on how to improve the efficiency of on-axis directional emission from PCWs. Some of them [3–5] still followed the original model in Moreno’s report [1], whereas some others proposed totally new strategies such as introducing two point defects near the PCW termination [6], covering a self-collimation photonic crystal (PC) at the termination of PCW [7], adding a coupled-resonator optical waveguide near the input or termination of PCW [9]. Besides the study of on-axis directional emission, there are also some recent works [12, 13] concerning the realization of off-axis directional emission from PCWs. For instance, Caglayan *et al*. [12] suggested a method to obtain off-axis directional beam by adding an asymmetrical gratinglike layer on the modified termination of PCW.

The aforementioned works play an important role in solving the problem of diffraction limit in normal PCW. However, there still remain two critical issues. The first one is whether the efficiency of directional emission in each previous work has reached the maximum. For case of Moreno’s model, it is well known that, the performance of directional emission mainly depends on the coupling efficiency between PCW and surface modes [5, 14], and such coupling efficiency is sensitive to the surface parameters such as radius, refractive index [3], and lattice period [13] of surface cylinders. In many previous works [1, 3–5, 12], the lattice period of cylinders in the modified termination or added surface was fixed at the same or twice as the period of bulk PC. So it is possible to further enhance the coupling efficiency if one can optimize the period of surface cylinders and other parameters simultaneously. The other issue is whether one can exactly control the beaming angle of directional emission. Though different methods have been proposed to realize off-axis directional emission [12, 13], yet no concrete formulation has been reported to exactly control the beaming angle.

Even if one can provide such formulation, it is still hard to solve the efficiency problem in issue one at the same time, because the coupling efficiency and beaming angle are both influenced by many surface parameters while a summative formulation to describe the relationships between them seems hardly feasible. So it is an arduous work if one still uses a step-by-step thought to analyze each parameter for solving the above issues. Instead, a reverse thought, without considering the aforementioned relationships, is employed in this work to solve these problems. Provided a surface added behind the PCW, a detector is placed at arbitrary location in the output field to characterize the transmission from waveguide [3]. Then, a design of highly-efficient directional emission with pre-defined beaming angle can be obtained by adopting an optimization algorithm to search the maximum power value of detector. In particular, we choose the Genetic Algorithm (GA) [15, 16] as optimization algorithm. As a global optimization algorithm, GA has the following features: robustness, capacity to search without gradient information, independence of the initial guess, and ability to operate on discrete and continuous parameters simultaneously. In general, GA is particularly useful when dealing with complex problems, in which the physical details are relatively unclear. Many previous works have demonstrated such effectiveness of GA in designing kinds of photonic structures [17–22]. Hence, it is very suitable to do the optimization task in our reverse thought.

## 2. Model and methods

#### 2.1 Description of the to-be designed system

As photonic band-gap (PBG) system, the two-dimensional (2D) PC considered here is composed of a square lattice of circular rods with a refractive index of *n*=3.4 (e.g., InGaAsP-InP at a wavelength of 1.55*µ*m) in vacuum. The radius of the dielectric rods is *r*=0.18*a*, where *a* is the lattice constant of bulk PC. In this paper, we only investigate the case of TE polarization (electric field pointing along the rod axis). For TE polarization, the PBG of the PC structure is calculated by the plane wave expansion (PWE) method and can be found in the frequency range between *ω*=0.30×2*πc/a* and *ω*=0.44×2*πc/a* [1, 3], where c is the light velocity in vacuum. As shown in Fig. 1, the directional emission system for either on-axis or off-axis case is based on a PCW, formed by removing a row of rods along the *x*-direction from the above mentioned 2D PC. The dimension of such waveguide is 10*a* in *x*-direction and 31*a* in *y*-direction.

To verify the reverse thought as previously said in the introduction, behind the termination of PCW, we add a gratinglike surface constructed by a column of cylinders made by the same dielectric as PCW. Compared with the modified terminations [1, 3–5] or added surfaces [10, 11] in previous works, the gratinglike surface [12, 13] has a simpler configuration, which can bring more convenience in further fabrication. It is expected to achieve highly-efficient directional emission by optimizing all possible parameters in gratinglike surface. For on-axis case (Fig. 1(a)), we realize that the optimum system must correspond to a symmetric structure, which is obvious from the symmetry of the problem. We divide the whole gratinglike surface into two symmetrical sub surfaces along the axis of y=0. The to-be optimized parameters are *r*
_{1}, *r*
_{2}, *r*
_{3}, and *r*
_{4}. Here, *r*
_{1} and *r*
_{2} represent the diameter and lattice period of dielectric rods in each sub surface, whereas *r*
_{3} and *r*
_{4} represent the distance of each sub surface to the axis of *y*=0 and the termination of PCW. For such an on-axis directional emission system, a Gaussian continuous wave source is launched at the point (*x*=0, *y*=0) with frequency *ω*
_{0}. The width of the Gaussian wave is equal to the wavelength of incident wave *λ*. It excites waves that propagate along the PCW and the waves are then emitted at the added surface. To characterize the directional emission, we place a detector, centered at (*x*=*x _{D}*,

*y*=0), with cross-sectional length of 6

*a*in the output field. The model for off-axis directional emission is a little more complicated. Considering the non-symmetry of the problem, we divide the whole gratinglike surface into two different sub surfaces. As suggested by Caglayan’s work [12], arranging different lattice periods for two sub surfaces is a very effective technique to control the beaming angle. On the other hand, we think setting different sizes and locations for surface cylinders in both sides is helpful to enhance the efficiency of off-axis directional emission. So the to-be optimized variables for off-axis case are added to 8 as indicated in Fig. 1(b). The location of

*x*

_{0}=9.18

*a*represents the outlet of PCW and

*θ*can represent the off-axis degree of detector centered at (

*x*,

_{D}*y*), as well as the pre-defined beaming angle as one need. In the following sections, we will apply the GA to optimize the above proposed structures for achieving highly-efficient directional emission.

_{D}#### 2.2 Genetic algorithm

GA belongs to a family of stochastic search algorithms called evolutionary computation, the optimization process of which is similar to Charles Darwin’s theory of evolution [15]. To clarify the implementation process of GA in our study, a flowchart is presented in Fig. 2. As mentioned in section 2.1, there are several variables need to be optimized in the directional emission models. To begin with the GA, these variables are encoded by binary strings and directly stored in a vector, called a chromosome. The initial population P0 is created by randomly yielding Nc chromosomes. Each chromosome therefore compactly represents a gratinglike surface to be simulated. For example, the i-th chromosome to represent the gratinglike surface in Fig. 1(b) can be expressed by the following *L*-length vector

Here, *C _{k}*=0 or 1, (

*k*=1, 2, 3,⋯,

*L*) can represent a binary bit;

_{1}

_{2}

_{3}

_{8}

*Sr*,

*Sr*,

*Sr*,⋯,

*Sr*are eight binary strings representing the information of

*r*

_{1},

*r*

_{2},

*r*

_{3},⋯, r. Specially, in this paper, we set the search space for the above eight variables as [0,

*a*], [

*a*, 2

*a*], [0, 2

*a*], [0, 2

*a*], [0,

*a*], [

*a*, 2

*a*], [0, 2

*a*], and [0, 2

*a*], respectively. Besides, considering the requirement in practical fabrication, we assume the precision for all variables is 0.01

*a*. So the corresponding length of string for

*Sr*

_{1},

*Sr*

_{2},

*Sr*

_{3},⋯,

*Sr*

_{8}is 7, 8, 8, 8, 7, 8, 8, and 8, respectively. After encoding, each chromosome will be further decoded to determine the values of all variables. We utilize the following equation as decoding scheme

Here, ‘BTD’ is a function to translate the binary string to decimal value.

After decoding, the fitness value *F _{i}* of each chromosome will be evaluated by the fitness function. According to Ref. [3], the directed power of detector, normalized to the input power, obtained as time average by the finite-difference time-domain (FDTD) method [23], can characterize the efficiency of directional emission. So for the models in Fig. 1, we just use the following formulation as fitness function to determine the fitness value of each chromosome

Where

Here, **E**
_{S} and **H**
_{S} represent the electric and magnetic field intensity of light source, **E**
_{D} and **H**
_{D} represent the electric and magnetic field intensity of detector, *t*
_{1} and *t*
_{2} represent the beginning and stop time of detector, Δ*t* represents the time step of FDTD, and Δ*x* and Δ*y* represent the spatial interval of FDTD. In our study, to satisfy convergence of the FDTD technique, each periodic lattice in the computational domain is divided into 20×20 grid cells (Δ*x*=Δ*y*=0.05*a*), and the beginning and stop time of detector is set as 300 Δ*t* and 600 Δ*t*, respectively. So Ω(*ω*
_{0}, *t*) and *P*
_{D} (*ω*
_{0}) can represent the instantaneous and time-average normalized power of detector.

After fitness evaluation, the GA then works using the following genetic operators to evolve the population:

**Selection**. Selection creates a new generation by allotting more positions in the new population to those chromosomes with favorable fitness values, and eliminating those with poor values. Specially, we use the fitness-proportionate selection to choose chromosomes for further mating. In this selection scheme, a chromosome is selected out of *N _{c}* chromosomes with a probability

*P*that is proportional to its fitness

_{i}*F*, as shown in equation (5)

_{i}**Crossover**. Crossover is the basic operator for producing new chromosomes in the GA. Like its counterpart in nature, crossover produces new individuals that have some parts of both parent’s genetic material. Here, we adopt the simplest crossover operator that is single-point crossover in our study. In this procedure, a pair of parent chromosomes *S*⃗_{parent,1} and *S*⃗_{parent,2} are selected and mated to produce two child chromosomes *S*⃗_{child,1} and *S*⃗_{child,2} by taking a random convex combination of the parent vectors, as described by

Here, *h* is a random integer position from 1 to *L*-1 to represent the crossover point.

**Mutation**. In natural evolution, mutation is a random process where one allele of a gene is replaced by another to produce a new genetic structure. Similarly, in the GA, the role of mutation is often seen as providing a guarantee that the probability of searching any given string will never be zero and acting as a safety net to recover good genetic material that may be lost through the action of selection and crossover [16]. There are many kinds of mutation operators for binary GA. In particular, we use the simplest one, single-point mutation in our simulations. For an original chromosome *S*⃗_{orig} of length *L*, we select a random integer position, *q*, from 1 to *L* as the mutation point, and perform a logical ‘NOT’ operation on the mutation point to produce the mutated vector, *S*⃗_{mut}, as represented by

Now, the initial population has completed a generation of evolution. This process is repeated until the algorithm terminates, typically after a particularly “fit” solution is found, or more generally, after a pre-defined number of generations. The best chromosome in the final generation will be decoded to provide an optimum design of gratinglike surface.

## 3. Optimization results and discussion

In this section, we will present several GA designs for on-axis and off-axis directional emission. We initially assumed the frequency of the Gaussian continuous wave is equal to 0.35×2*πc/a*. Then three GA projects were executed for different to-be optimized structures, i.e., (Project 1) on-axis directional emission structure with the detector placed at *x _{D}*=50

*a*, (Project 2) off-axis directional emission structure with the detector placed at

*x*=50

_{D}*a*,

*θ*=10°, and (Project 3) off-axis directional emission structure with the detector located at

*x*=50

_{D}*a*,

*θ*=20°. In each GA project, we used population size of 200 chromosomes for each generation, and allowed the program to run for a total of 80 generations. The statistic information on searching process for each GA project is shown in Fig. 3, where the best fitness as a function of generation is indicated for each run of the GA. Our GA performed as expected, and we got a general increase of fitness as the algorithm progressed. The maximum fitness was obtained at 61, 56, and 58 generations for three projects, respectively. This steady performance of GA is due to its features of robustness, capacity to search without gradient information, and independence of the initial guess.

Specially, the optimum gratinglike surface for each project can be represented by: (Project 1) *r*
_{1}=0.62*a*, *r*
_{2}=1.61*a*, *r*
_{3}=1.47*a*, and *r*
_{4}=1.16*a*; (Project 2) *r*
_{1}=0.77*a*, *r*2=1.87*a*, *r*
_{3}=0.84*a*, *r*
_{4}=1.01*a*, *r*
_{5}=0.69*a*, *r*
_{6}=1.06*a*, *r*
_{7}=0.69*a*, and *r*
_{8}=1.26*a*; (Project 3) *r*
_{1}=0.36*a*, *r*
_{2}=1.49*a*, *r*
_{3}=0.45*a*, *r*
_{4}=1.31*a*, *r*
_{5}=0.66*a*, *r*
_{6}=1.07*a*, *r*
_{7}=1.38*a*, and *r*
_{8}=1.18*a*. Fig. 4 shows the spectrum of transmission versus wavelength when the normal incidence TE-polarization Gaussian wave is launched into the above three optimum PC structures. From Fig. 4, it can be seen that a maximum detector power of up to 0.970, 0.839, and 0.883 can be achieved when the frequency is 0.35×2*πc/a*, 0.351×2*πc/a*, and 0.349×2*πc/a*, respectively, for three optimum structures. We define an acceptable bandwidth as within which the detector power of the directional emission is larger than 80% of the maximum power. With this definition, from Fig. 4, it is found that the acceptable bandwidth for three optimum structures is applicable in [0.345×2*πc/a*, 0.352×2*πc/a*], [0.345×2*πc/a*, 0.355×2*πc/a*], and [0.345×2*πc/a*, 0.353×2*πc/a*], respectively. To simplify the discussion, in the following sections, we will restrict our attention to the frequencies of 0.35×2*πc/a*, 0.351×2*πc/a*, and 0.349×2*πc/a* for three optimum structures. It is worth noting that the transmission properties are similar to all frequencies in the acceptable bandwidth.

In order to verify the emission properties of our designed systems, the FDTD method with perfectly matched layers, which are absorption boundaries, is employed for calculating the propagation of the electromagnetic wave. The spatial interval of FDTD is still 0.05*a*. Fig. 5 displays the steady-state amplitude distribution of the electric fields for several structures. In particular, Fig. 5(a) shows the distribution for the case of normal PCW without any added surface. The frequency used here is 0.35×2*πc/a*. It is found that light exiting from the waveguide undergoes a strong angular spread. Fig. 5(b) represents the case of PCW with added gratinglike surface for on-axis design, i.e., the optimum structure in Project 1. Compared with Fig. 5(a), we can find that the angular divergence of the output beam is significantly narrowed and a directional emitting beam in the output field can be observed. Though the power detector is predefined at location of *x _{D}*=50

*a*in GA project, the output light can propagate a longer distance. The effective directional emitting beam terminates at around

*x*=70

_{D}*a*. Besides, there exist some lateral lobes of radiation in the output field, while the intensity of them is obviously weaker than the main directional emitting beam. Similar off-axis directional emitting beams and weak lateral lobes of radiation can also be observed in Figs. 5(c) and 5(d), which are corresponding to the optimum structures in Project 2 and Project 3, respectively.

To further investigate the directional emission properties for all structures in Fig. 5. We plot the angular distribution of the output intensity as shown in Fig. 6, which is obtained by placing a number of angular detectors in the output field (seen in Fig. 5(c)). The distance between the angular detectors and the outlet of PCW is *R*. Specially, we set *R*=40*a* for all calculations. As can be seen from Fig. 6(a), the output intensity exhibits a broad angular distribution for normal PCW. But when a gratinglike surface is added behind the termination of PCW, obviously, a nice directional distribution of the intensity is observed in Fig. 6(b). The energy is mainly concentrated in a very small range of about -3° to 3° with central degree of 0°, whereas the intensity of the lateral lobes of radiation is quite smaller than the main lobe. Similar phenomena can be observed from Figs. 6(c) and 6(d), where the energy is mainly concentrated in range of around 5° to 15° with central degree of 10° and in range of about 10° to 30° with central degree of 20°, respectively. These behaviors are in good agreement with our original design objectives.

Most importantly, to determine the efficiency of directional emission in our designs, the time-averaged transmission coefficients of power (normalized to the incident power) along the x-direction are calculated by placing the detector with width of 6*a* at different positions, e.g., the location of transmission detector in Fig. 5(b) varies from *x _{D}*=20

*a*to

*x*=70

_{D}*a*with step of 5

*a*. The calculated results are given in Fig. 7, from where we can see the transmission coefficient for case of on-axis design maintains at around 0.8 when

*x*varies from 20

*a*to 70

*a*, which indicates a high efficiency of such on-axis directional emission. At the same time, the transmission coefficients for two off-axis designs seem in good superposition and both decrease slowly as

*x*increases. For instance, for the case of off-axis degree

*θ*=10°, it can be clearly found that the transmission coefficient is above 0.7 within the range of 20

*a*≤

*x*≤50

*a*. When

*x*reaches a larger value of

*x*=70

*a*, the transmission coefficient only reduce by about 18%. These behaviors prove such off-axis emitting beam also own a high efficiency in the propagation direction.

At last, we want to give some discussions on the physical mechanism of highly-efficient directional emission in our designs. In Moreno’s model, the directional emitting beam arises from a constructive interference of light emitted from surface modes in corrugated termination. Based on this concept, the following authors [3–5] focused their efforts on enhancing the coupling efficiency between PCW and surface modes, as well as increasing the radiation from surface modes to obtain highly-efficient directional emission. The above explanation can not be directly used for our designs because our model is just a simple gratinglike surface without introducing any difference between odd and even cylinders. As indicated in previous works [1, 14], the surface mode in such kind of one-line added surface can not directly couple to the radiating mode of the free space. So it is believed that the directional beams in our designs are due to the interference of light emitted from both surface modes and the outlet of waveguide. At the same time, the interference is sensitive to many factors such as the coupling efficiency between PCW and surface modes, the intensity of radiation from surface modes, and the phase matching between all light. Since GA is a global optimization algorithm, which can automatically optimize all the above factors on directional emission, we can also obtain highly-efficient directional beams from a simple gratinglike surface.

## 4. Conclusion

In summary, we have proposed a reverse thought to realize highly-efficient directional emission from 2D PCW. In this thought, a gratinglike surface was added behind the termination of PCW, without considering the effects of surface parameters on directional emission, we have realized both on-axis and off-axis highly-efficient directional beaming effects by using a GA to maximize the power value of detector placed in the output field at different pre-defined locations. The enhanced transmission and beaming of light property could have some interesting applications in enhancing the light coupling efficiency from PCW to conventional waveguides or conventional optical fibers. In addition, for the fabrication reason, there are only four or eight surface parameters to be optimized in our model. It is possible to further enhance the efficiency of directional emission by optimizing the radius and location of each cylinder in the gratinglike surface individually, or trying other added surface, e.g., populating pixels of a rectangular [17] or round [18] grid (binary 0 or 1) behind the waveguide either with air or with the dielectric.

## Acknowledgements

This work was supported by doctorial creative fund of Nanjing University of Science and Technology. Besides, we would like to acknowledge the support of project “Plasmonics and Photonic Crystals” from Nanyang Technological University.

## References and links

**1. **E. Moreno, F. J. García-Vidal, and L. Martín-Moreno, “Enhanced transmission and beaming of light via photonic crystal surface modes,” Phys. Rev. B **69**, 121402 (2004).
[CrossRef]

**2. **P. Kramper, M. Agio, C. M. Soukoulis, A. Birner, F. Müller, R. B. Wehrspohn, U. Gösele, and V. Sandoghdar, “Highly directional emission from photonic crystal waveguides of subwavelength Width,” Phys. Rev. Lett. **92**, 113903 (2004).
[CrossRef] [PubMed]

**3. **S. K. Morrison and Y. S. Kivshar, “Engineering of directional emission from photonic-crystal waveguides,” Appl. Phys. Lett. **86**, 081110 (2005).
[CrossRef]

**4. **W. R. Frei, D. A. Tortorelli, and H. T. Johnson, “Topology optimization of a photonic crystal waveguide termination to maximize directional emission,” Appl. Phys. Lett. **86**, 111114 (2005).
[CrossRef]

**5. **D. Gan, Y. Qi, X. Yang, J. Ma, J. Cui, C. Wang, and X. Luo, “Improved directional emission by resonant defect cavity modes in photonic crystal waveguide with corrugated surface,” Appl. Phys. B **93**, 849–852 (2008).
[CrossRef]

**6. **C.-C. Chen, T. Pertsch, R. Iliew, F. Lederer, and A. Tünnermann, “Directional emission from photonic crystal waveguides,” Opt. Express **14**, 2423–2428 (2006).
[CrossRef] [PubMed]

**7. **D. H. Tang, L. X. Chen, and W. Q. Ding, “Efficient beaming from photonic crystal waveguides via self-collimation effect,” Appl. Phys. Lett. **89**, 131120 (2006).
[CrossRef]

**8. **Y. L. Zhang, Y. Zhang, and B. J. Li, “Highly-efficient directional emission from photonic crystal waveguides for coupling of freely propagated terahertz waves into Si slab waveguides,” Opt. Express **15**, 9281–9286 (2007).
[CrossRef] [PubMed]

**9. **Z. H. Zhu, W. M. Ye, J. R. Ji, X. D. Yuan, and C. Zen, “Enhanced transmission and directional emission via coupled-resonator optical waveguides,” Appl. Phys. B **86**, 327–321 (2007).
[CrossRef]

**10. **Z. F. Li, A. Koray, and O. Ekmel, “Highly directional emission from photonic crystals with a wide bandwidth,” Appl. Phys. Lett. **91**, 121105 (2007).
[CrossRef]

**11. **Q. Wang, Y. P. Cui, C. C. Yan, L. L. Zhang, and J. Y. Zhang, “Highly efficient directional emission using a coupled multi-channel structure to a photonic crystal waveguide with surface modification,” J. Phys. D: Appl. Phys. **41**, 105110 (2008).
[CrossRef]

**12. **H. Caglayan, I. Bulu, and E. Ozbay, “Off-axis directional beaming via photonic crystal surface modes,” Appl. Phys. Lett. **92**, 092114 (2008).
[CrossRef]

**13. **H. B. Chen, X. S. Chen, J. Wang, and W. Lu, “Tunable beam direction and transmission of light using photonic crystal waveguide,” Physica B **403**, 4301–4304 (2008).
[CrossRef]

**14. **W. Śmigaj, “Model of light collimation by photonic crystal surface modes,” Phys. Rev. B **75**, 205430 (2007).
[CrossRef]

**15. **J. Holland, *Adaptation in Natural and Artificial Systems* (University of Michigan Press, Ann Arbor, MI, 1975).

**16. **D. E. Goldberg, *Genetic Algorithms in Search, Optimization and Machine Learning* (Addison-Wesley, New York, 1989).

**17. **L. F. Shen, Z. Ye, and S. L. He, “Design of two-dimensional photonic crystals with large absolute band gaps using a genetic algorithm,” Phys. Rev. B **68**, 0351091 (2003).
[CrossRef]

**18. **L. Sanchis, A. Håkansson, D. López-Zanón, J. Bravo-Abad, and J. Sánchez-Dehesa, “Integrated optical devices design by genetic algorithm,” Appl. Phys. Lett. **84**, 4460–4462 (2004).
[CrossRef]

**19. **E. Kerrinckx, L. Bigot, M. Douay, and Y. Quiquempois, “Photonic crystal fiber design by means of a genetic algorithm,” Opt. Express **12**, 1990–1995 (2004).
[CrossRef] [PubMed]

**20. **J. Smajic, C. Hafner, and D. Erni, “Optimization of photonic crystal structures,” J. Opt. Soc. Am. A **21**, 2223–2232 (2004).
[CrossRef]

**21. **R. P. Drupp, J. A. Bossard, D. H. Werner, and T. S. Mayer, “Single-layer multiband infrared metallodielectric photonic crystals designed by genetic algorithm optimization,” Appl. Phys. Lett. **86**, 081102 (2005).
[CrossRef]

**22. **A. Husakov and J. Herrmann, “Chirped multilayer hollow waveguides with broadband transmission,” Opt. Express **17**, 3025–3035 (2009).
[CrossRef]

**23. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Artech House, Boston London, 2000).