The finite element method applied to the study of twodimensional photonic crystals.
Abstract
Calculations of the photonic band structure, transmission coefficients, and quality factors of various twodimensional, periodic and aperiodic, dielectric photonic crystals by using the finite element method (FEM) are reported. The fundamental equations governing the propagation of electromagnetic waves in inhomogeneous media are revisited together with the boundary conditions required for each of the performed calculations. A detailed account of the eigenvalue and harmonic propagation analysis of the electromagnetic problem is reported for several periodic and finitelength structures. It is found that this method reproduces quite well previous results for these lattices obtained with the standard plane wave method with regards to the eigenvalue analysis (photonic band structure calculations). However, in contrast with frequency methods, the finite element method easily allows one to study the timeharmonic propagation of electromagnetic fields and, thus, to calculate the transmission coefficient of finite clusters in a natural way. Moreover, the advantages of using this real space method for structures of arbitrary complexity are also discussed. In addition, point defect cluster quality factor calculations are reported by means of FEM and they are compared to the ones obtained with the FDTD and Harminv methods. As a result, FEM comes out as an effective, stable, robust, and rigorous tool to study light propagation and confinement in both periodic and aperiodic dielectric photonic crystals.
I 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 Joannopoulos1995 (); Sakoda2001 (); cefe (); istrate2006 (). 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) (Angel1, ; Angel2, ), 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 nonlinear 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 twodimensional heterostructures, such us threedimensional waveguides made of a twodimensional photonic crystal core sandwiched between two layers of substrate that confines light by simple refraction index matching. Therefore, investigating twodimensional 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 method (PW) Joannopoulos1995 (); Yablo1991 (); Johnson2001 () 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 mpbman () 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.
Interestingly, a method that has received very little attention in this field even though is has been known and extensively used in other areas of Physics and Engineering is the finite element (FE) method fem01 (); Lin2007 (); Frei2004 (); Sopaheluwakan2003 (); Pomar2004 (). Indeed, 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. 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, the demand of this resource is quite insensitive to the presence of defects, which could render this technique very advantageous for studying disordered lattices Sopaheluwakan2003 (). Therefore, at the present status of the photonic band gap materials field, it is important to assess whether this method could be useful to investigate the properties of photonic crystals from both the fundamental and applied points of view. This is precisely one of the aims of this paper: we report a comprehensive study of both the band structure and transmission coefficients of various twodimensional photonic crystals using the finite element method. The other main goal is to study the feasibility of this method to analyze resonant cavities in photonic crystals and their associated quality factors.
The structure of the paper is as follows: in the next section we revisit the equations that describe the propagation of electromagnetic fields in twodimensional inhomogeneous media and, in particular, in photonic crystals. Also, the boundary conditions required for each particular simulation are presented. Section III summarizes the results for topologies with the symmetry of the square and triangular lattices. Section IV examines the dispersion relations for crystals with point defects for both the square and triangular lattices. Then, defect mode tunability is briefly discussed. Section V 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 VII, we state the conclusions of this work.
Ii 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 chapter) is described by Maxwell’s equations. In particular, for twodimensional 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
(1) 
where is the th component of the electric field at position , is the inhomogeneous relative dielectric constant of the photonic crystal, and with the angular frequency of the incident electric field and the speed of light in free space. In writing down equation (1), it has been assumed that the photonic crystal is nonmagnetic () and nonconducting (). Once equation (1) is solved, the timeharmonic electric and magnetic fields are easily calculated as
(2)  
(3) 
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
(4) 
where is the th component of the magnetic field at position . The timeharmonic electric and magnetic fields are easily calculated once equation (4) is solved and they are given by
(5)  
(6) 
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,
(7) 
where 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,
(8) 
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, , the only effect on the EM field is a change of its phase,
(9) 
and
(10) 
for TE and TM polarization, respectively. In these expressions, is a vector of the photonic crystal lattice and 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
(11)  
(12) 
where and are the initial values of the electric and magnetic fields at the boundaries, respectively, and 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 nonreflecting.
Iii 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 nonmagnetic ^{I}^{I}IThis applies to the rest of photonic crystals considered in the present work. The dielectric constant of the rods has a value of . The ratio , where is the radius of the cylinders and the lattice parameter, was taken as . This well known structure was first investigated by McCall and coworkers 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 mesh elements and periodic boundary conditions given by Bloch’s theorem were implemented. For the MPB calculation, a resolution of () grid elements were used and the dielectric constant was average over grid points. The resulting photonic band structure is depicted in Fig. 1. As it is customary, the dimensionless quantity 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 component of the electric field, , were also calculated at the , , and points of the 1BZ and they are shown also in figure 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 nonphysical boundaries. The component of the electric field was set to and 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 ^{II}^{II}IISimulations 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.
Another important class of twodimensional 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 Joannopoulos1995 (). For this reason it is important to test whether the FE method produces correct results for this nonorthogonal 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 was used for the dielectric constant of the nonempty part of the photonic crystal. The radius of the cylinders is . 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 mesh elements and periodic boundary conditions were implemented across the boundaries. For the transmittance calculation in the direction, a cluster comprising unit cells (dielectric rods) was divided into mesh elements. For the transmittance calculation in the direction, a cluster comprising unit cells was divided into mesh elements. The same boundary conditions as for the square lattice calculations described above were used. The photonic band structure, transmittances in the and directions, and field patterns for TEpolarized 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 TEpolarized EM waves that can be clearly seen in the transmission spectra for light propagating in both the and 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 direction that, however, are not present in the 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 direction.
Even though that is not the main purpose of this paper, 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 point for the ninth band for different discretizations of the square unit cell. In the FE case, meshes with , , , and elements were used, whereas for the calculations done with MPB, resolutions of () grid elements, () grid elements, () grid elements, and () grid elements were used. We took the eigenvalue calculated with MPB by using a resolution of and a mesh size of as the exact one. The result of this comparison can be seen in Fig. 3a. 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. 3b, 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.
Iv 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 Edo () or to reduce the lasing threshold of a certain laser emitter Sibilia (). Besides, they have been described as key elements for several applications in many areas of physics and engineering such as enhancing high directivity antennas Temelkuran (), designing low power consumption and highly tunable optical buffering devices Tanabe (); Baba (), or constructing new PCVCSEL lasers Song () as well as for biosensing applications Rajesh (), 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.
iv.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. McCall (). In the present context, we have studied this topology by FE method and compared it to the predictions of the PWE method Johnson2001 (). Figure 4 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 . The rod dielectric constant has been set to , as in the original McCall’s work. A supercell has been used for the FEM calculations reported in this section. It has been discretized into mesh elements. A element grid has been used for the corresponding MPB calculations. The real part of the zcomponent of the electric field for the defect mode located at was calculated at the point of the 1BZ. The main drawback of the supercell approach is the bandfolding effect: redundant bands of the unit cell are folded back times ( being the linear dimension of the supercell). This fact leads to larger computational times since the amount of eigenvalues to be solved grows . We have therefore solved eigenvalues for points in space for bands by both methods. Experimental results obtained by McCall and coworkers, MPB calculations, and FEM results fully agree as it is shown in Figs. 4 and 5.
Two modes make their way across the upper band gap region but only the defect mode depicted in Fig. 4a concentrates its energy around the missing rod.
The density of states (DOS) shows that perturbation of a single rod induces a series of sharp Dirac peaks centered at the frequencies where the defect state occurs. DOS calculations have been also computed via FEM by randomly sampling kpoints constrained to the irreducible portion of the IBZ for each lattice. In the long wavelength limit, this quantity clearly exhibits the linear behavior expected for propagation in an homogeneous twodimensional dielectric medium.
iv.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.
Fig. 4b shows the corresponding results for the photonic crystal parameters used by Smith et al. Smith () 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. 4a 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 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 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 Joannopoulos1995 (). The symmetry of such point defect modes is analyzed in detail in Kim ().
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 gapmidgap ratio can be easily tuned by adequately preselecting rod radii and dielectric contrast. For rods of dielectric contrast equal to in square and triangular lattices, one can maximize the gapmidgap ratio for TE polarization by just adjusting the rod radii to . Band diagrams and DOS calculations for radius rods for square and triangular lattice are reported in Fig. 5a and Fig. 5b. There, localized defect bands are strongly confined by a wider band gap. Moreover, these defect modes can be spatially translated into frequency by simply choosing the appropriate index contrast between lattice rods and point defect rod Joannopoulos1995 (). 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 Cox (); Shen (); Geremia (). These facts will be discussed elsewhere unpublished_im ().
V 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 offaxis direction because the crystal has still infinite length in this direction, but inasmuch as the inplane 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 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
(13)  
(14) 
where and are the initial values of the electric and magnetic fields at the boundaries, respectively, and is the propagation constant. A set of monochromatic plane waves is excited for each frequency in the simulation domain. This way, a point defect is introduced in the topologies characterized in previous sections. The resulting transmittance results are depicted in Fig. 6.
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 (Fig. 6). 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 and dielectric constant of , wherein its dimensions have been gradually modified, the energy stored in the defect mode centered around is exchanged with the upper band energy situated around . The quality factor has been calculated by means of the transmittance response for each crystal, since by definition
(15) 
where 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 Brent (), which combines bisection, regula falsi, and inverse quadratic interpolation methods for root finding. In contrast with the accepted wisdom, this method has proved to be highly reliable and more efficient than using FDTD, as we explain below in more detail. 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 7 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 midgap frequency and its FWHM gets narrower whenever increases.
Vi 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 and compare the Q factor calculations attained so far using FEM. For the present manuscript, the FDTD simulations were performed using MEEP FDTD (), a freely available implementation package of this method. The success of FDTD method is due to its flexibility and to its adaptability to irregular or aperiodic geommetries. More generally, the FDTD method can compute a large number of frequencies at once and even extract modes of the spectrum.
In FDTD space and time is 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 splitfield PML proposed by Berenger Berenger () has been used. Such an artificial region is needed so as to absorb outgoing spurious waves from the computational grid rather than reflecting them back into the photonic crystal. Moreover, 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 normalincident waves will not be the dominant ones and they will be mixed with reflected waves. MEEP implements Uniaxial Perfectly Matched Layers (UPML) UPML (); PML_problem () in order to absorb outgoing spurious waves from the computational. Along this artificial medium, the PML is expressed in terms of effective anisotropic and MEEP ().
For the present FDTD computation, 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. 8). 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 stored power divided by the loss power after one cycle:
(16) 
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 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 has been used. Sampling the electric field response with a period relative to the resonance frequency, results in a quasiperiodic step function that induces some uncertainties in the Q factor determination if (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 (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) harminv () 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 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.
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.

Square lattice r=a  Square lattice r=a  
Cluster size 
FEM  FDTD  Harminv  FEM  FDTD  Harminv 
5 
72  88  108  107  
7  376  347  722  723  
9  1800  1765  4310  4308  
11  8350  7674  24500  24559  

Triangular lattice r=a  Triangular lattice r=a  
5  97  103  134  126  
7  358  316  1280  1166  
9  1240  1209  11100  11238  
11  5270  5650  90200  83892 
Vii Conclusions
The present manuscript reports a comprehensive study of the photonic properties of several twodimensional 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 wellknown plane wave and FDTD methods, but it has many advantages not present in the others. In contrast with frequencyspace 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 acurate 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 an effective 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 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 an effective, stable, robust, and rigorous tool 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 metallodielectric photonic crystals.
Acknowledgements.
We would like to thank the Basque Government for financial support under the SAIOTEK 2012 programme (ref. SIGMA) and the Grupos Consolidados 2006–2012 programme (Grant No. IT33107).References
 (1) J. J. Joannopoulos, R. D. Meade, 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, Advanced Materials 15, 16791704, (2003).
 (4) E. Istrate and E. H. Sargent, Photonic crystal heterostructures, Reviews of Modern Physics78, 455481, (2006).
 (5) Angel J. GarciaAdeva, Band gap atlas for photonic crystals having the symmetry of the kagome and pyrochlore lattices, New Journal of Physics 8, (2006).
 (6) Angel J. GarciaAdeva, Band structure of photonic crystals with the symmetry of a pyrochlore lattice, Physical Review B. emph73, (2006).
 (7) E. Yablonovitch, T. J. Gmitter, and K. M. Leung, Photonic band structure: The facecenteredcubic case employing nonspherical atoms, Physical Review Letters 67, 22952299, (1991).
 (8) S. G. Johnson and J. D. Joannopoulos, Blockiterative frequencydomain methods for Maxwell’s equations in a planewave basis, Optics Express 8, no. 3, 173190 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX83173. See also http://abinitio.mit.edu/wiki/index.php/MIT_Photonic_Bands.
 (9) MPB online manual, http://abinitio.mit.edu/wiki/index.php/MPB_manual.
 (10) J. Jin, The finite element method in electromagnetism, (Wiley–IEEE press, New York, 2002).
 (11) M. C. Lin and R. F. Jao, Finite element analysis of photon density of states for twodimensional photonic crystals with inplane light propagation, Optics Express 15, 207 (2007).
 (12) W. R. Frei and H. T. Johnson, Finiteelement analysis of disorder effects in photonic crystals, Physical Review B 70, 1651161–11 (2004).
 (13) A. Sopaheluwakan, Defect states and defect modes in 1D photonic crystals (MSc Thesis, University of Twente, 2003).
 (14) J. L. GarciaPomar and M. NietoVesperinas, Transmission study of prisms and slabs of lossy negative index media, Optics Express 12, 2081–2095, (2004).
 (15) S.L. McCall, P.M. Platzman, R.Dalichaouch, David Smith and S.Schultz, Microwave propagation in Two Dimensional Dielectric Lattices, Physical Review Letters (67), 2017â2020, (1991).
 (16) B. Temelkuran, Mehmet Bayindir, E. Ozbay, R. Biswas, M. M. Sigalas, G. Tuttle and K. M. Ho, Photonic crystalbased resonant antenna with a very high directivity, Journal of Applied Physics 87, 603605, (2000).
 (17) Takasumi Tanabe, Masaya Notomi, Eiichi Kuramochi, Akihiko Shinya and Hideaki Taniyama, Trapping and delaying photons for one nanosecond in an ultrasmall highQ photoniccrystal nanocavity, Nature Photonics 1, 4952, (2007).
 (18) Toshihiko Baba, Takashi Kawasaki, Hirokazu Sasaki, Jun Adachi and Daisuke Mori, Large delaybandwidth product and tuning of slow light pulse in photonic crystal coupled waveguide, Optics Express 16, 92459253, (2008).
 (19) DaeSung Song, SeHeon Kim, HongGyu Park, ChangKyu Kim and YongHee Lee, Singlefundamentalmode photoniccrystal verticalcavity surfaceemitting lasers, Applied Physics Letters 80, 39013903, (2002).
 (20) Rajesh V.Nair and R Vijaya Progress, Photonic crystal sensors:An overview, Progress in Quantum Electronics 34, 89134, (2010).
 (21) C.Sibilia, T.M.Benson, M.Marciniak and T.Szoplik, Photonic crystals: Physics and technology, Springer, (2008).
 (22) 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, Journal of the Optical Society of America B 10, 314321, (1993).
 (23) Edo Waks and Jelena Vuckovic, Coupled mode theory for photonic crystal cavitywaveguide interaction, Optics Express 13, 50645073 (2005).
 (24) SeHeon Kim and YongHee Lee, Symmetry Relations of TwoDimensional Photonic Crystal Cavity Modes, IEEE Journal of Quantum Electronics 39, 10811085, (2003).
 (25) Steven J. Cox and David C. Dobson, Maximizing Band Gaps In TwoDimensional Photonic Crystals, SIAM Journal of Appllied Mathematics 59, 2108â2120, (1999).
 (26) L.F. Shen, Z. Ye, and S. He, Design of twodimensional photonic crystals with large absolute band gaps using a genetic algorithm, Physical Review B 68, 0351095, (2003).
 (27) J.M. Geremia, J. Williams, and H. Mabuchi, Inverseproblem approach to designing photonic crystals for cavity QED experiments, Physical Review E 66, 06660612, (2002).
 (28) Ardavan F. Oskooi, David Roundy, Mihai Ibanescu, Peter Bermel, J. D. Joannopoulos and Steven G. Johnson, MEEP: A flexible freesoftware package for electromagnetic simulations by the FDTD method, Computer Physics Communications 181, 687â702, (2010).
 (29) V. A. Mandelshtam and H. S. Taylor, Harmonic inversion of time signals, Journal of Chemical Physics 107, 67566769, (1997). Erratum, ibid. 109, 4128 (1998).
 (30) Allen Taflove and Susan C. Hagness, Computational Electrodynamics: The FiniteDifference TimeDomain Method, Artech, (2000)
 (31) R.P. Brent, Algorithms for Minimization without Derivatives, PrenticeHall, (1973).
 (32) Ardavan F. Oskooi, Lei Zhang, Yehuda Avniel and Steven G. Johnson, The failure of perfectly matched layers and towards their redemption by adiabatic absorbers, Optics Express 16, 1137617, (2008).
 (33) Z. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee A perfectly matched anisotropic absorber for use as an absorbing boundary condition,IEEE Transactions on Antennas and Propagation 43, 14601463, (1995).
 (34) J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114, 185200, (1994).
 (35) Imanol Andonegui and Angel J. GarciaAdeva, (unpublished).