Open Access
Add to BibSonomyAbstract
A critical assessment of the finite element (FE) method for studying two-dimensional dielectric photonic crystals is made. Photonic band structures, transmission coefficients, and quality factors of various two-dimensional, periodic and aperiodic, dielectric photonic crystals are calculated by using the FE (real-space) method and the plane wave expansion or the finite difference time domain (FDTD) methods and a comparison is established between those results. It is found that, contrarily to popular belief, the FE method (FEM) not only reproduces extremely well the results obtained with the standard plane wave method with regards to the eigenvalue analysis (photonic band structure and density of states calculations) but it also allows to study very easily the time-harmonic propagation of electromagnetic fields in finite clusters of arbitrary complexity and, thus, to calculate their transmission coefficients in a simple way. Moreover, the advantages of using this real space method in the context of point defect cluster quality factor calculations are also stressed by comparing the results obtained with this method with those obtained with the FDTD one. As a result of this study, FEM comes out as an stable, robust, rigorous, and reliable tool to study light propagation and confinement in both periodic and aperiodic dielectric photonic crystals and clusters.
© 2013 Optical Society of America
1. Introduction
Photonic crystals have generated a surge of interest in the last decades because they offer the possibility to control the propagation of light to an unprecedented level [1–4]. In its simplest form, a photonic crystal is an engineered inhomogeneous periodic structure made of two or more materials with very different dielectric constants. When an electromagnetic wave (EM) propagates in such a structure whose period is comparable to the wavelength of the wave, unexpected behaviors occur. Among the most interesting ones are the possibility of forming a complete photonic band gap (CPBG) [5,6], that is, a frequency range for which no photons having frequencies within that range can propagate through the photonic crystal (PC), to localize light by introducing several types of defects in the lattice, or enhancing certain non-linear phenomena due to small or anomalous group velocity effects. Unfortunately, all these nice features come at a price: the length scales required in order to fabricate a photonic crystal appropriate for operation at telecommunication frequencies are below the micron, so that ingenious innovations were required in order to actually fabricate such structures. Furthermore, fabricating devices based on these lattices is even more challenging from a technological point of view, especially in three dimensions. This has resulted in a lot of effort being devoted to investigating photonic crystal devices based on two-dimensional heterostructures, such us three-dimensional waveguides made of a two-dimensional photonic crystal core sandwiched between two layers of substrate that confines light by simple refraction index matching. Therefore, investigating two-dimensional photonic crystals is not a mere academic exercise but an important task, both at the fundamental and applied levels.
From a theoretical perspective, one of the most important venues of research is to develop numerical methods to solve Maxwell’s equations of an EM wave propagating in a PC that are reliable, fast, and capable of dealing with large systems as close as possible to the real ones employed for experiments. The reason why numerical methods are so important in this field is rooted on the fact that the predictions of Maxwell’s equations for these systems are in excellent qualitative and quantitative agreement with experiments and they are able to describe the ample phenomenology exhibited by these systems. The development of these numerical methods has undergone an evolution that mimics to a certain extent the one underwent by the techniques used to calculate the electronic band structure of semiconductors. Thereby, while early research efforts focused on ideal infinite periodic photonic crystals, for which techniques on the frequency space are ideally suited (such as the plane wave expansion (PWE) method [1, 7, 8] or algorithms based on the tight binding method), it was soon realized that the real interest of these systems is disordered photonic crystals. Frequency methods can be used to address the calculation of certain quantities (such as the photonic band structure or density of states) for some types of disorder, mainly point or line defects, by performing supercell calculations [9] but they are not very useful for studying disorder that is not localized (such as random displacements of the dielectric constituents of the photonic crystal or compositional disorder), which are of great practical interest, as these types of disorder are usually associated to the fabrication techniques themselves. Furthermore, frequency methods, in principle, are only applicable to infinite systems whereas real systems are always finite. For these reasons, people have started to pay attention during the last years to real space methods. These methods differ from the previous ones in that they work with a finite cluster which does not need to be periodic.
In the particular case of the real-space FE method, there are a number of reasons that seem to suggest that it could be applied well to the study of the propagation of EM waves in inhomogeneous media [10]. Just to cite some advantages, this method allows one to study geometries of arbitrary complexity, it can deal with frequency dependent dielectric functions (metallic inhomogeneous structures) in a natural way, discontinuities in the dielectric function are not especially detrimental for convergence of the method, and the quantities are already calculated in the stationary regime. The only shortcoming of the method is its extensive computer memory usage. However, the demand of this resource is quite insensitive to the presence of defects, which could render this technique very advantageous for studying disordered lattices [11]. Interestingly, not only has this method received relatively little attention in the area of photonic nanostructures (with some notable exceptions, see for example [12–21]) but, what is even more concerning, some authors in the field raise legitimate concerns about the usefulness of this method for dielectric photonic crystals and resonant cavities: for example, Villeneuve et al [22] state that FEM is not very reliable with regards to calculate important quantities such as quality factors of resonant cavities, whereas Oskooi et al [23] point out that the flexibility of the FE method comes at the price of software complexity that may render FEM not very useful for dielectric devices operating at infrared frequencies. Since there are commercial implementations of the FE method geared towards optics and photonics that offer a nice, albeit flexible, interface to the final user, it seems timely to undertake a comprehensive study of the capabilities of the FE simulation method with regards to dielectric photonic materials.
The goal of this manuscript is to assess whether the FE method is adequate as an all-purpose computational method for investigating the optical properties of both periodic and aperiodic photonic crystals and clusters and resonant cavities. The main results coming out from this study are twofold: firstly, it is shown that the FE method is a stable, robust, rigorous, and reliable tool to study light propagation in inhomogeneous dielectric photonic materials; secondly, it is found that, in some cases, FEM is better suited than other, more commonly used, computational methods to calculate extremely sensitive quantities. In particular, this is the case for high-Q factors: using real-space harmonic propagation transmittance computations —together with Brent’s method [24] for determining the width of the resonance— allows for a more efficient and reliable calculation than analyzing the Q-factor in these structures using a time domain transient propagation environment.
The structure of the paper is as follows: in the next section we revisit the equations that describe the propagation of electromagnetic fields in two-dimensional inhomogeneous media and, in particular, in photonic crystals. Also, the boundary conditions required for each particular simulation are presented. In Section 3, we briefly describe the three computational methods used throughout the manuscript. Section 4 summarizes the results for topologies with the symmetry of the square and triangular lattices. Section 5 examines the dispersion relations for crystals with point defects for both the square and triangular lattices. Then, defect mode tunability is briefly discussed. Section 6 addresses finite extent point defect PCs wherein issues such as the effect of supercell dimensionality in the formation of photonic localized states and, hence, the improvement of its quality factor are discussed. Finally, in section 8, we state the conclusions of this work.
2. The mathematics behind the propagation of EM waves in inhomogeneous media
The propagation of an electromagnetic wave through an inhomogeneous medium (such as the photonic crystals considered in this section) is described by Maxwell’s equations. In particular, for two-dimensional media, any electric or magnetic field can always be expressed as a linear combination of a transversal electric (TE) and a transversal magnetic fields (TM). In the first case, the electric field is perpendicular to the photonic crystal plane —whereas the magnetic field is constrained into this plane— and the sourceless Maxwell’s equations are reduced to a Hemholtz equation for the electric field given by
where Ez(r⃗) is the z-th component of the electric field at position r⃗, ∊r (r⃗) is the inhomogeneous relative dielectric constant of the photonic crystal, and k0 = ω/c with ω the angular frequency of the incident electric field and c the speed of light in free space. In writing down equation Eq. 1, it has been assumed that the photonic crystal is non-magnetic (μr = 1) and non-conducting (σ = 0). Once equation Eq. 1 is solved, the time-harmonic electric and magnetic fields are easily calculated as On the other hand, in the case of TM polarization, the magnetic field is perpendicular to the photonic crystal plane —whereas the electric field is constrained into this plane— and the sourceless Maxwell’s equations reduce to a Hemholz equation for the magnetic field given by where Hz(r⃗) is the z-th component of the magnetic field at position r⃗. The time-harmonic electric and magnetic fields are easily calculated once equation 4 is solved and they are given by An important aspect of any electromagnetic simulation is the use of appropriate boundary conditions at the interfaces. On the one hand, for the photonic band structure calculations reported below, it is necessary to implement boundary conditions that mimic an infinite simulation domain together with the periodicity of the photonic crystal lattice. An infinite medium is simulated by using perfect magnetic (PMC) or perfect electric conductor (PEC) boundary conditions that mirror the simulation domain. The former condition is used for TE polarization of the EM field and ensures that the component of the magnetic field tangent to the boundary is identically zero at the outer interface, that is, where n̂ is a unit vector perpendicular to the outer simulation domain surface at each point. The later is used for TM polarization of the EM field and ensures that the component of the electric field tangent to the boundary is identically zero at the outer interface, that is, The periodicity of the photonic crystal lattice is ensured by an adequate use of Bloch’s theorem at the boundaries of the photonic crystal unit cell. This theorem states that when the electric (magnetic) field propagates from one point in the PC to another one separated from the previous one by a lattice vector, R⃗, the only effect on the EM field is a change of its phase, and for TE and TM polarization, respectively. In these expressions, R⃗ is a vector of the photonic crystal lattice and k⃗ is the wavevector of the electromagnetic wave. On the other hand, for the transmittance calculations reported below, finite clusters in the direction parallel to the incident wave vector were used, whereas the clusters are of infinite extension in the perpendicular direction. In order to mimic such a material, PMC and PEC boundary conditions were used for the outer interfaces that limit the simulation domain in the direction perpendicular to the incident wave vector for TE and TM polarization, respectively. For the interfaces at which the EM wave enters and exits the cluster, the situation is a little more involved because it is necessary to avoid unphysical reflections due to the finite size of the cluster and, thus, perfectly matched layers were used at these interfaces. The equations that describe such boundaries are given by where E0z and H0z are the initial values of the electric and magnetic fields at the boundaries, respectively, and β = k0 is the propagation constant. The upper condition applies to TE polarization, whereas the lower one corresponds to TM polarization. If the electric field is an eigenmode of the boundary, the boundary is exactly non-reflecting.3. Computational methods
The central computational method in this work is the finite element (FE) method [10]. This method is routinely used in several branches of engineering to simulate physical phenomena. However, its use in fundamental research and, particularly, in optics, is not so extended [12–19, 22], even though, several features of this method seem to suggest that it could apply very well to the study of the propagation of EM waves in inhomogeneous media. In particular, this method allows one to study geometries of arbitrary complexity, it can deal with frequency dependent dielectric functions (metallic inhomogeneous structures) in a natural way, discontinuities in the dielectric function are not especially detrimental for convergence of the method, and the quantities are already calculated in the stationary regime. The only shortcoming of the method is its extensive computer memory usage. However, in contrast with other methods (most noticeably the plane wave one), the demand of this resource is quite insensitive to the presence of defects, which renders this technique even more advantageous for studying disordered lattices [11].
Even though there are freely available implementations of this method for dealing with electromagnetic problems [25–27], they are more focused on other aspects of electromagnetism, most noticeably magnetostatics, antennae design, electrostatics, and so on. For this reason, we decided to use a commercial implementation of the FE method, namely the COMSOL® multiphysics package [28]. The main reasons to choose this software are the comfortable CAD environment it provides for designing the system to be simulated, there is an specific application mode for simulating propagation of EM waves in an arbitrary medium (including both 2D and 3D systems) and its post-processing capabilities (it can be interfaced with MATLAB®).
At the core of any FE simulation is the method for generating a mesh, that is, a partition of the geometry into small units of known shape called mesh elements. In two dimensions, this program uses an unstructured mesh (triangular mesh elements) generator based on the Delaunay algorithm. Once a mesh is created, the dependent variables are approximated by a known function (shape functions) that can be described with a finite number of parameters called degrees of freedom. Inserting this approximation into the original equations generates a set of equations for the degrees of freedom that is then solved with an appropriate solver. In particular, for the present problem, the shape functions are second order Lagrange elements and the solver for the resulting linear problem is the so called UMFPACK solver, which is a very efficient direct solver for unsymmetric systems [29].
The PWE method is the most popular method for obtaining Bloch modes and eigenfrequencies in dielectric periodic systems such as photonic crystals. The PWE method poses Maxwell’s equations in a decoupled Hermitian form. Then, the fields are expanded in a Fourier series along the directions defined by the reciprocal lattice vectors, yielding a numerically solvable algebraic eigenvalue problem. The well-known MIT Photonic-Bands (MPB) [8] implementation of the PWE method has been used in this manuscript as the standard benchmark to assess the accuracy of the finite element method regarding the calculation of the photonic band structures and densities of states. Some of the nice features of this software package are that it is fully vectorial; it supports arbitrary, anisotropic dielectric structures and non-orthogonal unit cells; parallel computations are supported by using the MPI library; it provides a flexible user interface (extensible and scriptable); and it has been extensively tested by numerous researchers during the last years.
As mentioned above, the PWE method is not very well suited to deal with finite length or aperiodic systems and the band folding phenomenon makes it very cumbersome for analyzing defects in photonic crystal structures. In contrast, the FDTD method overcomes some of these limitations [30]. In fact, it is fair to say that this is the simulation method most widely used for studying the propagation of EM fields in photonic crystals and other inhomogeneous structures. FDTD performs a direct integration of the traveling wave through the PC using a dual discretization: both in time and space. This way, the input wave propagates through the structure by time stepping along the grid. In this work, the free implementation of the FDTD method known as MEEP [23] has been used in order to discuss the advantages and disadvantages of the FE method when dealing with defect cavities in photonic crystals and the corresponding estimation of their quality factors. As a bonus, MEEP includes an enhanced interpolation and subpixeling of the simulation domain and it supports nonlinear materials.
4. Simulation results for photonic crystals based on square and triangular lattices
The first structure analyzed in this work is a photonic crystal made of dielectric circles whose centers occupy the positions of a square lattice and is depicted in the inset of Fig. 1. The dielectric material was assumed to be linear, isotropic, and non-magnetic (this applies to the rest of photonic crystals considered in the present work). The dielectric constant of the rods has a value of 9. The ratio , where r is the radius of the cylinders and a the lattice parameter, was taken as 0.38. This well known structure was first investigated by McCall and coworkers [31] in order to compare the predictions of theory with experimental results with regards to the localization of light in strongly scattering media. In the present context, we have studied this topology in order to check how well the FE method fares when compared with the plane wave method that, as stated above, is the one commonly used to calculate photonic band structures. The first nine photonic bands were calculated for transversal electric (TE) polarization along the path that delimits the irreducible part of the 1st Brillouin Zone (1BZ). For the FE calculation, the square unit cell was divided in 3720 mesh elements and periodic boundary conditions given by Bloch’s theorem were implemented. For the MPB calculation, a resolution of 64 ×64 (=4096) grid elements were used and the dielectric constant was average over 9 grid points. The resulting photonic band structure is depicted in Fig. 1. As it is customary, the dimensionless quantity ωa/2πc = a/λ has been used to characterize the frequency of the incident EM wave, where ω is the frequency of the incident EM wave and λ the associated wavelength. The corresponding eigenmodes of the z-component of the electric field, Ez, were also calculated at the Γ, M, and K points of the 1BZ and they are shown also in Fig. 1. It is clear from that figure that the band structure calculated with the FE method faithfully reproduces the one calculated with MPB to its minimum details. There are three photonic gaps in the band structure of this lattice whose sizes coincide with the calculated ones with MPB. Also, the modes calculated with the FE method closely resemble those calculated with MPB up to a trivial symmetry operation or linear combination of degenerate modes. In addition, a quantity that can be readily calculated with the FEM method is the transmittance of a finite photonic crystal. The transmittance of the cluster is calculated using the usual approach of power integration along a domain boundary. The boundary conditions for the simulation domain were set as follows: on the input and output boundaries (left and right boundaries) a perfectly matched layer was used to avoid spurious reflections from non-physical boundaries. The z-component of the electric field was set to 1 and 0 at the initial time of the simulation on the input and output boundaries, respectively. On the other hand, perfectly magnetic conductor boundary conditions that set the tangential magnetic component of the magnetic field to zero were used for the upper and lower boundaries (simulations using different boundary conditions for the upper and lower boundaries were performed and no significant differences in the calculated transmission spectra were found) in order to mimic an infinite stripe in that direction.
Fig. 1 On the left: band structure calculated with the MPB (blue circles) and COMSOL® (red line) software packages for TE-polarized EM waves. On the right: Ez patterns in the unit cell calculated with MPB (upper row in each rectangle) and COMSOL® at the special symmetry points Γ, X, and M.
Another important class of two-dimensional photonic crystals is the one formed by those structures based on the triangular lattice because this is the symmetry most commonly used for practical applications [1]. For this reason it is important to test whether the FE method produces correct results for this non-orthogonal lattice. A photonic crystal made of dielectric cylinders whose centers occupy the sites of an hexagonal lattice of lattice constant a (see the insets in Fig. 2) was analyzed. A value of 12 was used for the dielectric constant of the non-empty part of the photonic crystal. The radius of the cylinders is 0.12a. The simulation setup was very similar to the one used for the square lattice: For the FE band calculation, the unit cell was divided into 1086 mesh elements and periodic boundary conditions were implemented across the boundaries. For the transmittance calculation in the ΓM direction, a cluster comprising 12 ×8 unit cells (dielectric rods) was divided into 60748 mesh elements. For the transmittance calculation in the ΓK direction, a cluster comprising 11 ×10 unit cells was divided into 67282 mesh elements. The same boundary conditions as for the square lattice calculations described above were used. The photonic band structure, transmittances in the ΓM and ΓK directions, and Ez field patterns for TE-polarized EM waves are reported in Fig. 2. The band structure calculated with the FE method faithfully reproduces the one obtained with MPB. This structure possesses one gap between the first and second bands for TE-polarized EM waves that can be clearly seen in the transmission spectra for light propagating in both the ΓM and ΓK directions as wide opaque regions. There are also some partial gaps, as the ones occurring between the second and third bands and the fourth and fifth bands in the ΓM direction that, however, are not present in the ΓK direction. Also, there are some opaque regions not related to gaps on the band structure but rather to the existence of uncoupled modes, such as the one associated to the 5th band, that is uncoupled in the ΓK direction.
Fig. 2 Left panel: Transmittance for TE-polarized waves propagating in the ΓM direction of the hexagonal photonic crystal cluster depicted in the inset. Central panel: Comparison between the band structure calculated with COMSOL® (red line) and MPB (blue circles) for TE-polarized EM waves along the boundary of the irreducible part of the 1BZ. The inset shows the unit cell used for the calculations. Right panel: Transmittance for TE-polarized waves propagating in the ΓK direction of the hexagonal photonic crystal cluster depicted in the inset. Inside the green rectangle: Comparison between the Ez field patterns at the M point calculated with MPB (on the left) and FEM (on the right) for the first eight bands (the band index increases from bottom to top). A portion of the photonic crystal that contains 3 ×3 unit cells is displayed in order to show the hexagonal symmetry of the modes.
It is illustrative to compare to some extent the accuracy and speed of the FE calculations with the ones performed by using the PW method. To estimate the relative accuracy, we computed the percent error in the eigenvalue calculated at the X point for the ninth band for different discretizations of the square unit cell. In the FE case, meshes with 254, 414, 928, and 1502 elements were used, whereas for the calculations done with MPB, resolutions of 16 ×16 (=256) grid elements, 32 ×32 (= 1024) grid elements, 48 ×48 (= 2304) grid elements, and 64 ×64 (= 4096) grid elements were used. We took the eigenvalue calculated with MPB by using a resolution of 256 ×256 and a mesh size of 25 as the exact one. The result of this comparison can be seen in Fig. 3(a). It is noticeable that the FE method gets a better accuracy with coarser discretizations of the lattice than the PW does. This is due to the use of second order Lagrange elements. However, the differences between both methods should be negligible for most applications. Of course, this better accuracy comes at a price and the simulation times for the FE calculations are longer by a noticeable factor than those for the MPB ones, as can be seen in Fig. 3(b), where we have depicted the evolution of the simulation runtime with the number of mesh elements and grid elements for the FE and MPB calculations, respectively.
Fig. 3 (a) Evolution of the relative error with the number of elements in which the lattice is discretized. (b) Evolution of the simulation run time with the number of elements in which the lattice is discretized.
5. Defect states in 2D photonic crystal superlattices
In complete analogy with semiconductors, the physics of disordered photonic crystals is even richer and their potential applications even more promising. Among these, photonic crystals in which defects are placed in a controlled way forming point or line defects are especially important for telecommunication applications. Their importance is rooted on the fact that these types of defects can lead to the formation of localized states inside the photonic gap that allows one to localize light around the defects.
Point defect photonic crystal configurations have been widely studied as promising candidates for the enhancement of strong coupling between the resonance cavity and quantum dots light emitters [32] or to reduce the lasing threshold of a certain laser emitter [33]. Besides, they have been described as key elements for several applications in many areas of physics and engineering such as enhancing high directivity antennas [34], designing low power consumption and highly tunable optical buffering devices [35,36], or constructing new PC-VCSEL lasers [37] as well as for biosensing applications [38], just to cite some examples.
In the simulations reported in this section, periodicity of the lattice has been intentionally broken and hence a defect is introduced into the otherwise perfectly ordered dielectric distribution, giving rise to localized states within the band gap region. When a single rod is removed from the dielectric lattice, light bounces back and forth in the disordered area, trapped by the surrounding band gap, whilst no other leakage mechanism is present. The resulting photonic crystal can still be classified on the basis of unit cell calculations, as long as the fundamental periodic cell hosts sufficient rods around the imperfection to make sure that the defect state does not overlap with neighboring replicas of itself. Furthermore, in these single rod perturbed cells, rotational symmetry remains invariant and thus comparative studies with previous perfect square and triangular lattice based calculations can still be performed at the edge of the 1BZ.
5.1. A square lattice of dielectric rods with a point defect: the McCall’s experiment
The experimental evidence of spatial mode localization in an ordered dielectric lattice of rods in an air background was first given by McCall et al [31]. In the present context, we have studied this topology by FE method and compared it to the predictions of the PWE method [8]. Figure 4(a) depicts the resulting dispersion diagram calculated for the experimental setup raised by McCall, where a single rod has been suppressed amidst the square lattice of dielectric rods, for which the normalized rod radius is 0.38a. The rod dielectric constant has been set to 9, as in the original McCall’s work. A 5 ×5 supercell has been used for the FEM calculations reported in this section. It has been discretized into 95903 mesh elements. A 3096 element grid has been used for the corresponding MPB calculations. The real part of the z-component of the electric field for the defect mode located at a/λ =0.4686 was calculated at the Γ point of the 1BZ.
Fig. 4 (a) On the left: Band structure of a 5 ×5 square lattice periodic supercell calculated with MPB (blue circles) and FEM (red lines), wherein a defect state has been excited around = 0.466 ±0.002. On the right: density of photonic states (DOS), calculated by FEM. The band gap region is clearly distinguishable and a weakly localized defect mode merging from the upper band can be seen. In the inset, Ez patterns for a 5 × 5 supercell calculated by both methods are shown for the M symmetry point. (b) On the left: band structure of a 5×5 triangular lattice calculated with MPB (blue circles) and FEM (red line). A single defect state merges at = 0.285 ±0.003. On the right: density of photonic states, calculated by FEM. In the inset, Ez patterns for the Γ symmetry point. An hexagonal supercell has been used for the FEM calculation.
The main drawback of the supercell approach is the band-folding effect: redundant bands of the unit cell are folded back N times (N being the linear dimension of the supercell). This fact leads to larger computational times since the amount of eigenvalues to be solved grows ∼ N. We have therefore solved eigenvalues for 211 points in K-space for 100 bands by both methods. Experimental results obtained by McCall and coworkers, MPB calculations, and FEM results fully agree, as shown in Fig. 4(a).
Density of states (DOS) calculations were also computed via FEM by randomly sampling k-points constrained to the irreducible portion of the 1BZ for each lattice and they are also reported in Fig. 4(a). In the long wavelength limit, this quantity clearly exhibits the linear behavior expected for propagation in an homogeneous two-dimensional dielectric medium. The band gaps are also clearly visible as the frequency regions of zero DOS.
5.2. A triangular lattice of dielectric rods with a point defect
An analogous situation can be achieved by removing a single rod in a photonic crystal based on the triangular lattice. Figure 4(b) shows the corresponding results for the photonic crystal parameters used by Smith et al. [39] in their investigation of the defect mode structures in the square and triangular lattices. Rod radii and dielectric constant parameters of the triangular lattice supercell are the same as in the orthogonal case of Fig. 4(a) and, once again, FEM calculation and results obtained by means of PWE method coincide. The simulation setup was similar to the previously described one for the square lattice, but this time discretization of the supercell has been set to 71376 mesh elements. Furthermore an hexagonal supercell has been used for the FEM calculations, as its shape matches the reciprocal lattice of a periodical triangular arrangement. In both cases, square and triangular lattice point defect supercells, the electric field z-component is well localized around the defect neighborhood and rapidly decays to small amplitudes as one moves further from the defect. Also, the eigenfunction around the lattice irregularity clearly shows an inherent symmetry. Indeed, both lattices present a monopole pattern with a single nodal plane through each dielectric rod [1]. The symmetry of such point defect modes is analyzed in detail in [40].
The confinement of the electric field amplitude in the cavity is directly related to the band gap strength. Stronger mode localization has been observed for wider gap geometries. This is a highly desirable effect since the gap-mid-gap ratio can be easily tuned by adequately preselecting rod radii and dielectric contrast. For rods of dielectric contrast equal to 9 in square and triangular lattices, one can maximize the gap-mid-gap ratio for TE polarization by just adjusting the rod radii to 0.17a. Band diagrams and DOS calculations for 0.17a radius rods for square and triangular lattice are reported in Fig. 5(a) and Fig. 5(b), respectively. There, localized defect bands are strongly confined by a wider band gap and the DOS diagrams show a sharp peak centered at the resonance frequency. In contrast, the structures depicted in Figs. 4(a) and 5(b) possess a less localized resonant mode near the upper edge of the photonic band gap and thus, they are not as clearly seen in the corresponding DOS diagram. Moreover, these defect modes can be spatially translated into frequency by simply choosing the appropriate index contrast between lattice rods and point defect rod [1]. With regards to this mode tunability some important properties of photonic crystal cavities can be conveniently enhanced. Indeed, incremental design procedures have successfully been demonstrated by means of algorithmic techniques rather than using intuitive trial and error design methods, where the defect dielectric function is tuned according to design requirements [41–43]. These facts will be discussed elsewhere [44].
Fig. 5 (a) On the left: band structure of a 5 ×5 square lattice periodic supercell calculated with MPB (blue circles) and FEM (red lines). The rod radii has been set to 0.17a, which ensures a better localization of the unique defect mode. On the right: density of photonic states, calculated by FEM. The defect mode clearly shows up in the gap region, since the defect state is strongly confined. In the inset, the Ez patterns for a 5 ×5 supercell calculated with both methods are shown for the M symmetry point. All these patterns match the ones presented in Fig. 4(a) but the modal volume decreases significantly for 0.17a rod lattices. (b) Triangular lattice made of 0.17a rods where the central rod has been removed. The band gap increases and so it does the defect mode localization around the point defect.
6. Defect states in finite 2d photonic crystal clusters
All the structures studied so far had an infinite extension in all space directions. Nevertheless, an engineered photonic crystal must be finite sized. This fact limits defect mode confinement factor as stored energy will be finally radiated outside the cavity. This circumstance gives rise to a quantity of very much practical importance: the quality factor of a photonic crystal resonator, Q, which is just the ratio between the stored energy in the photonic crystal and the radiated energy per cycle. In the present context it is unnecessary to consider the radiated energy in the off-axis direction because the crystal has still infinite length in this direction, but inasmuch as the in-plane propagation will be limited to N periods of dielectric rods, electromagnetic fields will still be leaked outwards. Quantitatively, these energy decay mechanisms are uncorrelated to each other, and so they can be studied separately. On the other hand, in the present calculations no losses, absorptions, nor imperfections have been taken into account and, consequently, the quality factor is only limited by the inherent energy leakage from the cavity to the radiative modes. Transmittance diagrams have been obtained using the FEM but, in contrast with the transmittances reported above for ideal lattices, in this section a finite length design is regarded. Therefore, PMC and PEC boundary conditions were substituted by PML interfaces in the direction perpendicular to the incident wave vector as well as at the interfaces at which the EM wave enters and exits the structure. Accordingly, transparent influx conditions were imposed to the outer boundaries of the finite cluster and a TE polarized source propagating in the Γ-X direction was considered. Spurious reflections are thus prevented from being injected back into the simulation domain by locating these artificial perfectly matched layers all around the simulation domain boundaries. The equations that describe such boundaries are given by
where E0z and H0z are the initial values of the electric and magnetic fields at the boundaries, respectively, and k⃗0 is the propagation constant. A set of monochromatic plane waves was excited for each frequency in the simulation domain. The resulting transmittance results are depicted in Figs. 6 and 7 for square and triangular arrangements of dielectric rods, respectively. As can be seen in those figures, when the size of the cluster is enlarged by adding subsequent periods to the lattice, the transmittance spectra clearly reveals a sharp peak centered at the defect mode frequency. When one considers a larger cluster, the defect mode gets narrower, i.e. it converges to the previously discussed periodic case. Additionally, the defect mode frequency is shifted due to the overall dielectric percentage change in the finite cluster which is indeed very convenient since it allows one to tune the resonance frequency according to dielectric index changes. For the particular case of a square arrangement of dielectric rods with radius 0.38a and dielectric constant of 9, wherein its dimensions have been gradually modified, the energy stored in the defect mode centered around = 0.4686 is exchanged with the upper band energy situated around = 0.49. The quality factor has been calculated by means of the transmittance response for each crystal, since by definition where fc represents the central frequency of the defect mode and FWHM the difference between the two values in the transmittance function at which the transmittance is equal to the half of defect mode maximum amplitude. For each quality factor calculation, the FWHM has been accurately determined by means of the Brent’s algorithm [24], which combines bisection, regula falsi, and inverse quadratic interpolation for root finding. The method of computing Q-factors by directly obtaining data from the transmittance spectrum has been reapeatedly critiziced [22, 23], mainly due to the low efficiency and instabilities presented by time domain methods for performing iterative transmittance data computations. In contrast, FEM is very efficient and reliable in calculating transmittance spectra, which yields an efficient way of computing accurate Q values from Eq. 15. This is illustrated in more detail in the next section, where we compare the Q values obtained by using FEM with those obtained by using FDTD and harmonic inversion of time signals methods.Fig. 6 (a) Transmittance curves for different cluster sizes in a rectangular arrangement of dielectric rods of radius 0.38a with a central defect. (b) Analogous to (a) for radius 0.17a. The insets show a detailed view of the defect mode curves indicated by the red arrow. Logarithmic scale has been used for the vertical axis.
Fig. 7 Analogous to figure 6 but for triangular arrangements of dielectric rods with a central defect. The radii of the rods are (a) 0.38a and (b) 0.17a.
As the amount of dielectric rods surrounding the defect increases, the rate of the energy loss within the cavity relative to the energy confined in it decays exponentially. Figure 8 depicts the quality factor increment for these square and triangular clusters. According to that figure, the defect mode bandwidth decreases exponentially with the addition of new periods due to the strengthening of the band gap effect, which is the main contribution to the enhancement of the quality factor. The onset of the leakage mechanism transforms the localized modes transmittance spectra into distorted Lorentzian peaks. These peaks progressively tend to a Lorentzian curve when the defect mode gets closer to the mid-gap frequency and its FWHM gets narrower whenever N increases.
Fig. 8 Top: quality factor estimation as a function of the Lorentzian shaped resonance FWHM for a number of different sized square clusters of dielectric rods with radii (a) 0.38a and (b) 0.17a. The FWHM value is given in dimensionless normalized units of a/λ. Logarithmic scale has been used for both axes. Down: analogous computation results obtained for triangular arrangements of rods with radii (c) 0.38a and (d) 0.17a.
7. Time domain approach: FDTD
As FDTD is the most popular method for the theoretical description of light propagation in these systems, we used FDTD in order to assess the reliability of the Q factors calculated using FEM. For the present manuscript, the FDTD simulations were performed using MEEP [23], a freely available implementation of this method. The success of the FDTD method is due to its flexibility and to its adaptability to irregular or aperiodic geometries. More generally, the FDTD method can compute a large number of frequencies and extract the modes of the spectrum at once, which renders it the method of choice for calculating high-Q factors in PC topologies. However, computing Q-factors using Eq. 15 and FDTD can be cumbersome. In FDTD one needs to compute the transmittance spectrum at the end of the computational cell, somewhere before the PML, running the simulation until the fields evolve and decay significantly (|Ez|2 or |Hz|2 need to decay by at least a factor of ). Then, it is necessary to normalize this flux by running the simulation again without the resonant cavity structure. Furthermore, the time needed to compute the direct transmission simulation grows with the value of the quality factor of the resonance cavity since, according to Fourier’s theorem, the number of time steps scales as the inverse of the frequency resolution one needs to attain [23]. Therefore, computing Q-factors by means of the transmittance spectrum obtained from FDTD is not a fast, nor reliable, approach.
Even if a time domain transient propagation analysis is used to calculate the Q-factors (rather than the transmittance data directly), there are still other complications associated to using FDTD. Indeed, in FDTD space and time are divided into a finite rectangular grid and then fields are evolved in time using discrete steps. However, FDTD has serious limitations when the computational domain is finite whilst FEM convergence time is rather insensitive to this fact. In FDTD the computational domain must be terminated with some boundary conditions as it is the case in the FEM and, in order to simulate open boundaries in finite clusters, a wave absorbing mechanism, such as the split-field PML proposed by Berenger [45] has been used. Such an artificial region is needed in order to absorb outgoing spurious waves from the computational grid rather than reflecting them back into the photonic crystal. Unfortunately, quality factor calculations are very sensitive to the size of the computational grid and, thus, if the injected pulse experiences spurious reflections from the domain boundaries, radiated normal-incident waves will not be the dominant ones and they will be mixed with reflected waves.
As briefly mentioned above, the FDTD method easily allows for the possibility of either using broadband or monochromatic sources, being the former the most effficient way to excite the resonance. FEM, on the other hand, can only deal with monochromatic sources in the stationary harmonic regime. Therefore, when comparing the FEM results of the previous section with the ones calculated with FDTD, there is the question of whether using broadband or monochromatic excitation for the FDTD calculations. To address this point it is important to take into account that performing iterative (i.e., single frequency) transmittance data computations by using time domain methods is very inefficient and subject to instabilities [22, 23]. Therefore, comparing FEM to FDTD using this type of excitation for both methods, would set the FDTD in a disadvantageous situation from the onset. For this reason and since we are interested in comparing the best possible Q-factors results obtained by both methods, a broadband source has been used for FDTD simulations, whereas a monochromatic excitation was used for the FEM ones, as explained above. In particular, for the present FDTD computations, a point Gaussian source has been placed inside the cavity. This excitation must be short enough (broad bandwidth) to excite the defect mode for each cluster. When this Gaussian source is switched on, the field grows and after some time this source is extinguished. Subsequently, a resonance effect occurs and the electromagnetic fields bounce back and forth for a limited amount of optical periods. Meanwhile, the energy trapped around the defect exhibits an exponential time decay (see Fig. 9). At this point, if one takes into account the slowly varying component only, i.e., the envelope of the electric field norm, the quality factor is determined by the ratio of the power stored at a given time t divided by the power lost after one period T
Fig. 9 (a) evolution of the electromagnetic field inside the defect rod for a 0.38a arrangement and for different cluster dimensions. A point source is excited in the cavity but after some periods the source is extinguished and although the electromagnetic field still remains in the cavity it experiences an exponential decay. (b) analogous computation results obtained for triangular defect cluster.
Regarding the PML implementation used for the present FDTD simulations, we have used Uniaxial Perfectly Matched Layers (UPML) [46, 47] implemented in MEEP in order to absorb outgoing spurious waves from the computational domain. Along this artificial medium, the PML is expressed in terms of effective anisotropic ∊ and ν[23]. It is important to notice that, in this type of calculations, the proximity of the UPML to the cluster significantly determines the decaying factor of the existing resonance. In the FDTD calculations reported in Table 1 an UPML of thickness 2a surrounds each structure. Besides, the distance between the UPML layer and the edge of the computational grid must be adequately tuned in order to avoid an unphysical field decay. According to this, a distance of 0.8 ≤ d ≤ 1.1 has been used. Sampling the electric field response with a period relative to the resonance frequency, results in a quasi-periodic step function that induces some uncertainties in the Q factor determination if Eq. 16 is directly applied. These fluctuations could be due to the abruptly broken translational symmetry of the clusters. However, the marked decaying behavior of the electric field data can be easily filtered and interpolated. Therefore, by iteratively applying Eq. 16 to the interpolated sample field, an average Q factor and the corresponding theoretical error is obtained. Finally, in order to compare with previous Q factor calculations obtained via FEM transmittance results and FDTD decaying field value relations, an harmonic inversion of time signals ( Harminv) [48] implementation has been used. This software package uses a filter diagonalization method (FDM) to find the deconvolution of sinusoidal functions near a given frequency interval. The simulation setup is the same as the one used in the standard FDTD calculations: a broad Gaussian point source located in the defect excites TE modes in the cavity and PML layers surrounding the cluster. In particular, a = 0.02 pulse width Gaussian source centered at the resonant frequency of each structure has been excited. When the source vanishes, Harminv performs the signal processing of the fields in the cavity. This way, it identifies the frequencies and decay rates of the excited resonant modes. This method is a fast and accurate way to determine the Q factors but still requires waiting until the fields have evolved and decayed to a certain small value. Computation time scales also with the cluster size, but this fact can often be overcome if symmetries are used to reduce the size of the computational cell. Besides, choosing a broad Gaussian source can yield to the appearance of spurious solutions but, at the same time, the source must be broad enough to excite the resonant frequency.
Table 1. Quality factor calculations for different cluster sizes in square and triangular lattice rod-type PC single-defect cavities using FEM, FDTD, and harmonic inversion methods.
It is important to stress the remarkable agreement between FEM results, FDTD and Harminv calculations of the Q factor. Table 1 summarizes the results obtained using the three methods described in this report. Two dimensional FEM calculation of Q factors offers a very satisfactory agreement with FDTD and Harminv results whilst providing substantial gain in solution robustness and efficiency. Indeed, FDTD is a powerful tool to calculate resonant frequencies and quality factors for complex cavity structures. However, it is very inefficient because one must discard many simulation cycles before reaching the stationary regime. Also, there are many heuristic factors that enter into the FDTD simulation process, such as the thickness of the UPML, the excitation source location and both its central frequency and width, that make the results from this method not as reliable as one would want. Moreover, when computing resonant modes in time domain especial care must be taken to avoid choosing an excitation source nearly orthogonal to the resonant mode because FDTD is likely to miss it, otherwise. In this regard, FEM can be seen as an efficient, reliable and more rigorous alternative to FDTD for the analysis of quality factors and resonant modes in complex dielectric material.
8. Conclusions
The present manuscript reports a comprehensive study of the photonic properties of several two-dimensional photonic crystals and finite clusters by using the finite element method. The main result coming out from these calculations is that the FEM allows one to reproduce the results obtained with the well-known plane wave and FDTD methods, but it has many advantages not present in the others. In contrast with frequency-space based approaches, FEM also deals in a natural way with finite clusters and aperiodic structures of arbitrary complexity. To prove this, we have calculated the band structure of periodic photonic crystals based on the square and triangular lattices. It is found that the band structures calculated in this way are almost indistinguishable of those calculated with the well known MPB package. Also, the modes calculated with FEM closely resemble those calculated with the PWE method. It is noticeable that using a coarser discretization FEM results are more accurate than the ones given by PWE.
Moreover, the transmission coefficients of a number of finite clusters of the aforementioned lattices were calculated along different directions in reciprocal space and for both TE and TM polarizations. The features (such as the position and width of the photonic gaps) of these transmittances agree quite well with the band structures and once again the accuracy of the band and mode calculations is very good when compared with the PWE method. Therefore, these results demonstrate that the FEM method can be a very useful general purpose method for investigating photonic crystals. This fact is further stressed by the results obtained for PCs containing a single defect, where DOS and dispersion diagrams have been calculated. There, experimental results obtained by McCall and coworkers, MPB calculations, and FEM results fully agree. As seen in these cavities, the confinement of the electric field amplitude strongly depends on the band gap strength of the underlying structure and, thus, wider gap geometries support higher mode localization. In addition, former calculations have been reproduced for finite length point defect clusters. In this context, the localization of light around the defect region has been quantified by accurately determining the quality factor using FEM, FDTD, and Harminv procedures for different cluster dimensions.
FEM is proven to be a reliable and stable tool for point defect cluster quality factor calculation, wherein for each quality factor calculation, the FWHM has been accurately determined by means of Brent’s algorithm. This technique permits to determine the essential information needed for Q factor calculation in a speedy and computationally effective way. In addition, the leakage mechanism of PC cavities transforms the transmittance spectra into an almost Lorentzian peak and, therefore, with few transmittance calculations one can obtain the entire point defect transmittance response. It is noteworthy, with regards to the the point defect PCs addressed above, the fact that the three methods accurately reproduce the Q factor for each topology. However, FDTD based methods have some serious drawbacks. On the one hand, the width and the proximity of the absorbing layers significantly determine the decaying behavior of the trapped fields inside the cavity, and so, the parameters of the PMLs must be carefully chosen for each structure. In the aforementioned point defect structures, the source must be placed inside the cavity, which in some cases can be an unrealistic situation. Furthermore, the bandwidth of the source must be set intuitively so as to excite only the defect mode. On the other hand, we believe that among these methods FEM is desirable for the calculation of Q factor due to its high numerical efficiency and stability because in FDTD, after the source is extinguished, one must wait an uncertain amount of time until the fields evolve and decay. In fact, FDTD gives quite accurate values for both the resonant frequency and the Q factor but, for higher Q values, the slope of the electromagnetic field inside the defect is nearly zero and hence, convergence time and and numerical errors increase drastically.
To summarize, the results presented in this manuscript demonstrate that the finite element method is a stable, robust, rigorous, and reliable tool that perfectly complements other, more extended, techniques in order to study light propagation and confinement in both periodic and aperiodic dielectric photonic crystals. Furthermore, we expect that these advantages can be extrapolated to systems in which the optical constants are frequency dependent, such as hybrid metallo-dielectric photonic crystals.
Acknowledgments
We would like to thank the Basque Government for financial support under the SAIOTEK 2012 ref. ( S-PE12UN043) programme. One of us, I.A., wants to thank the Vicerrectorado de Euskara y Plurilinguismo of the UPV/EHU for financial support under the PhD Fellowships 2011 programme.
References and links
1. J. J. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, New Jersey, 1995).
2. K. Sakoda, Optical Properties of Photonic Crystals (Springer, Berlin, 2001).
3. C. Lopez, “Material aspects of photonic crystals,” Adv. Mater. 15, 1679–1704, (2003). [CrossRef]
4. E. Istrate and E. H. Sargent, “Photonic crystal heterostructures,” Rev. Mod. Phys. 78, 455–481, (2006). [CrossRef]
5. A. J. Garcia-Adeva, “Band gap atlas for photonic crystals having the symmetry of the kagome and pyrochlore lattices,” New J. Phys. 8, 86/1–14 (2006). [CrossRef]
6. A. J. Garcia-Adeva, “Band structure of photonic crystals with the symmetry of a pyrochlore lattice,” Phys. Rev. B. 73, 0731071 (2006). [CrossRef]
7. E. Yablonovitch, T. J. Gmitter, and K. M. Leung, “Photonic band structure: the face-centered-cubic case employing nonspherical atoms,” Phys. Rev. Lett. 67, 2295–2299, (1991). [CrossRef] [PubMed]
8. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001). [CrossRef] [PubMed]
9. MPB on-line manual, http://ab-initio.mit.edu/wiki/index.php/MPB\_manual.
10. J. Jin, The Finite Element Method in Electromagnetism (Wiley–IEEE press, New York, 2002).
11. A. Sopaheluwakan, Defect States and Defect Modes in 1D Photonic Crystals (MSc Thesis, University of Twente, 2003).
12. M. -C. Lin and R. -F. Jao, “Finite element analysis of photon density of states for two-dimensional photonic crystals with in-plane light propagation,” Opt. Express 15, 207–218 (2007). [CrossRef] [PubMed]
13. W. R. Frei and H. T. Johnson, “Finite-element analysis of disorder effects in photonic crystals,” Phys. Rev. B 70, 1651161 (2004). [CrossRef]
14. J. L. Garcia-Pomar and M. Nieto-Vesperinas, “Transmission study of prisms and slabs of lossy negative index media,” Opt. Express 12, 2081–2095, (2004). [CrossRef] [PubMed]
15. M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O. Painter, “Optomechanical crystals,” Nature 462, 78–82, (2009). [CrossRef] [PubMed]
16. M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O. Painter, “A picogram- and nanometre-scale photonic-crystal optomechanical cavity,” Nature 459, 550–555, (2009). [CrossRef] [PubMed]
17. A. H. Safavi-Naeini, T. P. M. Alegre, M. Winger, and O. Painter, “Optomechanics in an ultrahigh-Q two-dimensional photonic crystal cavity,” Appl. Phys. Lett. 97, 1811061–3, (2010). [CrossRef]
18. E. Gavartin, R. Braive, I. Sagnes, O. Arcizet, A. Beveratos, T. J. Kippenberg, and I. Robert-Philip, “Optomechanical coupling in a two-dimensional photonic crystal defect cavity,” Phys. Rev. Lett. 106, 2039021 (2011). [CrossRef]
19. Huang-Ming Lee and Jong-Ching Wua, “Transmittance spectra in one-dimensional superconductor-dielectric photonic crystal,” J. Appl. Phys. 107, 09E1491–3, (2010).
20. V. F. Rodriguez-Esquerre, M. Koshiba, and H. E. Hernandez-Figueroa, “Finite-Element Analysis of Photonic Crystal Cavities: Time and Frequency Domains,” J. Lightwave Technol. 23, 1514–1521, (2005). [CrossRef]
21. J. K. Hwang, S. B. Hyun, H. Y. Ryu, and Y. H. Lee, “Resonant modes of two-dimensional photonic bandgap cavities determined by the finite-element method and by use of the anisotropic perfectly matched layer boundary condition,” J. Opt. Soc. Am. B 15, 2316–2324, (1998). [CrossRef]
22. P. R. Villeneuve, S. Fan, and J. D. Joannopoulos, “Microcavities in photonic crystals: mode symmetry, tunability, and coupling efficiency”, Phys. Rev. B 547837–7842, (1996). [CrossRef]
23. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: a flexible free-software package for electromagnetic simulations by the FDTD method,” Comp. Phys. Comm. 181, 687–702, (2010). [CrossRef]
24. R. P. Brent, Algorithms for Minimization without Derivatives (Courier Dover Publications, 1973).
25. Elmer – Finite Element Software for Multiphysical Problems, http://www.csc.fi/elmer/index.phtml.
26. The unofficial numerical electromagnetic code archives, http://www.si-list.org/swindex2.html.
27. The EMAP Finite Element Modeling Codes, http://www.emclab.umr.edu/emap.html.
28. Comsol multiphysics and Electromagnetics module, http://www.comsol.com.
29. T. A. Davis, UMFPACK 4.6: Unsymmetric MultiFrontal sparse LU factorization package, http://www.cise.ufl.edu/research/sparse/umfpack/.
30. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).
31. S. L. McCall, P. M. Platzman, R. Dalichaouch, D. Smith, and S. Schultz, “Microwave propagation in two dimensional dielectric lattices,” Phys. Rev. Lett. 67, 2017–2020, (1991). [CrossRef] [PubMed]
32. E. Waks and J. Vuckovic, “Coupled mode theory for photonic crystal cavity-waveguide interaction,” Opt. Express 13, 5064–5073 (2005). [CrossRef] [PubMed]
33. C. Sibilia, T. M. Benson, M. Marciniak, and T. Szoplik, Photonic Crystals: Physics and Technology (Springer, Milano, 2008). [CrossRef]
34. B. Temelkuran, M. Bayindir, E. Ozbay, R. Biswas, M. M. Sigalas, G. Tuttle, and K. M. Ho, “Photonic crystal-based resonant antenna with a very high directivity,” J. Appl. Phys. 87, 603–605, (2000). [CrossRef]
35. Takasumi Tanabe, Masaya Notomi, Eiichi Kuramochi, Akihiko Shinya, and Hideaki Taniyama, “Trapping and delaying photons for one nanosecond in an ultrasmall high-Q photonic-crystal nanocavity,” Nat. Photonics 1, 49–52, (2007). [CrossRef]
36. T. Baba, T. Kawasaki, H. Sasaki, J. Adachi, and D. Mori, “Large delay-bandwidth product and tuning of slow light pulse in photonic crystal coupled waveguide,” Opt. Express 16, 9245–9253, (2008). [CrossRef] [PubMed]
37. D. -S. Song, S. -H. Kim, H. -G. Park, C. -K. Kim, and Y. -H. Lee, “Single-fundamental-mode photonic-crystal vertical-cavity surface-emitting lasers,” Appl. Phys. Lett. 80, 3901–3903, (2002). [CrossRef]
38. R. V. Nair and R. Vijaya, “Photonic crystal sensors: an overview,” Prog. Quant. Electron. 34, 89–134, (2010). [CrossRef]
39. D. R. Smith, R. Dalichaouch, N. Kroll, S. Schultz, S. L. McCall, and P. M. Platzman, “Photonic band structure and defects in one and two dimensions,” J. Opt. Soc. Am. B 10, 314–321, (1993). [CrossRef]
40. S. -H. Kim and Y. -H. Lee, “Symmetry relations of two-dimensional photonic crystal cavity modes,” IEEE J. Quantum Electron. 39, 1081–1085, (2003). [CrossRef]
41. S. J. Cox and D. C. Dobson, “Maximizing band gaps in two-dimensional photonic crystals,” SIAM J. Appl. Math. 59, 2108–2120, (1999). [CrossRef]
42. L. F. 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, 035109 (2003). [CrossRef]
43. J. M. Geremia, J. Williams, and H. Mabuchi, “Inverse-problem approach to designing photonic crystals for cavity QED experiments,” Phys. Rev. E 66, 066606 (2002). [CrossRef]
44. I. Andonegui and A. J. Garcia-Adeva, (unpublished).
45. J. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phys. 114, 185–200, (1994). [CrossRef]
46. A. F. Oskooi, L. Zhang, Y. Avniel, and S. G. Johnson, “The failure of perfectly matched layers and towards their redemption by adiabatic absorbers,” Opt. Express 16, 113761 (2008). [CrossRef]
47. Z. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Ant. Prop. 43, 1460–1463, (1995). [CrossRef]
48. V. A. Mandelshtam and H. S. Taylor, “Harmonic inversion of time signals,” J. Chem. Phys. 107, 6756–6769, (1997). Erratum, ibid. 109, 4128 (1998). [CrossRef]
