## Abstract

We demonstrate the capabilities of principal component analysis (PCA) for studying the results of finite-difference time-domain (FDTD) algorithms in simulating photonic crystal microcavities. The spatial-temporal structures provided by PCA are related to the actual electric field vibrating inside the photonic microcavity. A detailed analysis of the results has made it possible to compute the phase maps for each mode of the arrangement at their respective resonant frequencies. The existence of standing wave behavior is revealed by this analysis. In spite of this, some numerical artifacts induced by FDTD algorithms have been clearly detailed through PCA analysis. The data we have analyzed are a given set of maps of the electric field recorded during the simulation.

©2004 Optical Society of America

## 1. Introduction

Since the late 1980s, after the publication of some seminal papers (for example, Ref. [1]), the number of experimental and numerical results on photonic crystal structures have increased. Among the wide variety of structures analyzed, photonic crystal microcavities [2] have become a new tool for enclosing radiation. Some nontrivial phenomena—the modification of spontaneous emission, the Kerr effect in photonic crystals [3, 4]—have been identified, and great effort has been devoted to investigating and clarifying the causes of these facts. Numerical simulations are essential tools for obtaining a deeper insight into the physics underlying photonic crystal behavior. One approach that is often used is to solve Maxwell equations through a finite-difference time-domain numerical scheme (FDTD) [5]. Time-domain computations allow for an unambiguous characterization of important properties of microcavities, such as band-gap features, resonant defect mode dynamics, and far-field pattern properties of crystal modes. The results have frequently been interpreted by use of fast Fourier transform analysis and filtering techniques [6]. A complete numerical characterization of the photonic crystal is often required, guiding the design process. Furthermore, other numerical tools and optimization algorithms can be incorporated into this process to improve and refine the final results [7].

Therefore, techniques for extracting useful information about mode characteristics should be available to photonic crystal researchers. In this paper, we propose applying a technique based on principal component analysis (PCA) to analyze the results obtained from FDTD simulations of two-dimensional photonic crystal microcavities. PCA has been successfully employed for analysis in image-processing systems and in noise characterization of cameras [8]. Here we illustrate how PCA methods produce a more accurate determination of photonic crystals properties.

Section 2 explains the PCA analysis. In Section 3 we focus on the photonic crystal microcavity that we have simulated. We have chosen a well-known system: a two-dimensional photonic crystal with a central defect [6]. The cavity is made of GaAs cylinders in air. It allows the presence of localized modes in the bandgap: the resonant modes of the system. They are excited selectively by placing a specific set of point sources in different configurations. The three subsections developed here are devoted to the analysis of specific features of the electric field distribution that remained unnoticed without the application of the PCA method. A spatial distribution of the mode’s phase is computed in Subsection 3.1 as a consequence of the grouping of couples of principal components with the same harmonic temporal evolution but shifted *π*/2 between them. Subsection 3.2 analyzes the meaning of another principal component that appears systematically after the PCA action. Its electric field distribution is related to the existence of a standing wave within the microcavity. Another interesting result derived from PCA is the capability to discern the presence of artifacts resulting from symmetry mismatch. This capability is demonstrated in Subsection 3.3. Finally, Section 4 summarizes the main findings of the paper.

## 2. Basis of the principal component analysis

PCA has demonstrated its utility for the classification and analysis of several sets of images. It is based on multivariate statistical techniques [9]. In our case, the input data are a set of *N* frames containing maps of a given electromagnetic component regularly spaced in time, {*F*
_{1},*F*
_{2},…,*F _{N}*}, where

*F*=

_{i}*F*(

*t*), with

_{i}*t*the instants of time where the frames are taken. Each frame contains

_{i}*M*points arranged as one-, two-, or three-dimensional maps. If the map is not one dimensional, it is reshaped as a one-dimensional vector by a procedure that preserves the information about the actual location of each spatial point. The statistical variables are the electromagnetic component at different time instants,

*F*. They configure an

_{i}*N*-dimensional variable, {

*F*}. The

_{i}*M*points define

*M*realizations of the previously defined statistical variable. The large number of statistical realizations usually available when analyzing irradiance maps or electromagnetic field maps has made possible a clearer interpretation of the obtained data. Actual data exhibit nonnegligible correlations among frames. PCA generates a collection of new frames, or variables, that are uncorrelated. They are the so-called principal components,

*Y*, where

_{α}*α*runs from 1 to

*N*. These new variables are obtained by a diagonalization procedure applied to the covariance matrix,

*S*. The covariance matrix is calculated between frames. Then,

*S*is a

*N*×

*N*matrix. The diagonalization procedure can be written as follows,

where *λ* denotes the eigenvalues of this equation that represent the contributions of the principal component with the same label, *α*; and *e _{α}* denotes the eigenvectors of the previous diagonalization. In addition, these eigenvectors are used to generate the principal components from the original data,

This last equation shows how the principal components can be written as a linear combination of the original frames. In some manner, the PCA method can be seen as a linear change of coordinates. Moreover, the previous equation can be inverted to obtain the original frames from the principal components. When doing such a reconstruction, we can also limit or select the principal components. The result will be a filtered version of the original frames. This is also known as principal component rectification.

The results of the PCA applied to *N* frames are *N* eigenvalues (*λ _{α}*),

*N*eigenvectors (

*e*), and

_{α}*N*principal components (

*Y*). Each principal component can be arranged as a map having

_{α}*M*points. We may define them as eigenframes or eigenimages. The eigenvectors are orthogonal and represent the temporal evolutions of the contribution of each eigenimage to the original data. When describing frames with a given power spectrum density for a temporal-stationary phenomena, PCA can be used to sample the spectrum and characterize it. The involved principal components have quasi-harmonic time dependence [10]. Finally, the eigenvalues, explain in decreasing order the contribution of their associated principal component to the total variance of the original data. It is also interesting to note that the capability of the PCA for sectioning the variance and quantify the individual contributions of the obtained principal components to the total variance of the original data, has made possible to reveal spatial-temporal structures hidden behind the main contributions to the variance of the data. To quantify this, we define the following parameter,

where *λ _{α}* denotes the eigenvalues obtained from the PCA, and

*w*is the relative contribution of the given principal component,

_{α}*α*, to the total variance of the data.

When we characterize noise within this multivariate technique, it is possible to associate an uncertainty to each eigenvalue [9]. This uncertainty is used to connect two consecutive eigenvalues when their respective uncertainties overlap [8, 11]. Then a process is defined as the frames retrieved when only the principal components associated with consecutively overlapped eigenvalues are used. The concept of process is used to ease the interpretation of the results obtained from the PCA. The classification into processes reduces and groups the meaningful principal components and provides an analytical tool that can be automatically applied and implemented.

FDTD data are generated by computational algorithms that have been scrutinized and refined over time to improve their output and make them more accurate and reliable. However, their implementation will be affected by intrinsic limitations and constraints derived from the finite and discretized realizations of the algorithm, both in time and in space. PCA has been proposed as a tool to identify and quantitatively estimate the presence of numerical artifacts and noise in the FDTD output. An important feature of FDTD data is that the “signal” is much larger than the “noise.” In this case the signal is the primary source of the variance of the data. This fact, along with the decreasing importance of ordering of the eigenvalues, will identify the actual electromagnetic field distributions with the first principal components. These “signal” principal components may also be grouped into several processes according to the previously established method for defining a process. These features will be used and illustrated in the following sections.

## 3. Photonic crystals analyzed by PCA

Our goal in this paper is to demonstrate the utility of PCA through a detailed study of a “benchmark” example. A complete description of the simulated structure can be found in Ref. [6]. For a better understanding we have followed the same choices and simulation parameters as those used in Ref. [6]. The analyzed photonic crystal can support five modes, two of which are degenerated. The resonance frequencies can be found in Ref. [6]. Simulations are two dimensional. The nonzero components of the electromagnetic fields are *E _{z}*,

*H*, and

_{x}*H*(this situation corresponds to a TMz mode). We are going to excite some modes of the cavity in order to study them with the PCA technique. By using this method, we will discriminate processes associated with different variance structures and we will ascribe them to relevant physical quantities of the modes. At this point, we should note this fact: PCA methods can be applied to any sequence of frames and will be efficient in the analysis of any other photonic crystals. The dimensional limitations are more related to the actual capabilities of the computing hardware. The PCA method is a blind tool that analyzes the structure of the variance of the data. Here the source excitation methods, the size and location of the computational grid, or the geometry and nature of the spatial structure will not compromise the application of PCA. However, all these constitutive parameters and properties are necessary for properly understanding and analyzing the PCA results.

_{y}In the following subsection we will illustrate some of the most relevant features that have been obtained from the results of the PCA. In Fig. 1 we have represented the FDTD results for the monopolar mode of the simulated geometry. The computational spatial grid is composed of 221 × 221 nodes, the Courant factor is *s* = 0.7068, and the spatial step is *δ* = 0.025*μ*m. We obtained the data by selecting the electric field maps every 10 temporal steps. This rate guarantees that the temporal evolution is correctly sampled. The number of frames in the sequence of data is 101. In this case, the starting point for recording the data is placed after 40.000 temporal steps. These results are in perfect accordance with the referenced results of Guo and Albin [6]. Figure 2 represents the eigenvalues in a semilog plot. At the same time, we have applied the classification method to group together those principal components with eigenvalues that are undistinguished. In this figure, we have also included the maps of some eigenimages (numbers 20,40,60, and 80) belonging to the same process. This process makes up 0.00328% of the total variance. The temporal evolutions of the first principal components are shown in Fig. 3. The quasi-harmonic dependencies are clearly visible. In this case, eigenvectors *e*
_{1} and *e*
_{2} oscillate as harmonic functions and are shifted *π*/2 one with respect to the other. This frequency is the same as the monopolar mode frequency. In representing *e*
_{1}(*t*) versus *e*
_{2}(*t*), the points are located on a circle, thus validating the existence of the quasi-monochromatic process. The analysis of the eigenimages is left for the following subsections.

#### 3.1. Phase map computation

The application of PCA for the analysis of noise figures made possible the definition of spatial-temporal processes that were obtained during analysis of the uncertainty and distinguishability of the eigenvalues obtained in the PCA application. These processes reduce the complexity of the variance structure by grouping together those structures that have eigenvalues with consecutive overlapping uncertainties. The uncertainties themselves were obtained with statistical techniques. On the other hand, the temporal evolution of the variance spatial structures, the eigenimages, are described by the associated eigenvectors. This fact was studied more in depth for finally relating the eigenvalues with sampled values of the power spectral density of a given sequence of frames [10]. In that case, the eigenvectors showed a clear harmonic evolution that allowed a temporal frequency interpretation. In the case analyzed in this paper, we have found some other connections between the temporal evolution of the eigenvectors. These findings have been interpreted within the scope of the electromagnetic computations performed by the FDTD methods. Then we have proposed grouping together those eigenimages associated with eigenvectors that have the same frequency in their quasi-harmonic temporal evolution. As shown below, this association is made in couples. Moreover, their temporal evolutions are shifted *π*/2 one with respect to the other. This is the only solution for the eigenvectors to remain orthogonal at the same frequency [10]. Therefore, we may propose the existence of “quasi-monochromatic processes” formed by the eigenvectors (and their associated eigenimages and eigenvalues) with the same temporal frequency.

A relevant feature obtained from the PCA for the photonic crystal study is that two principal components, *Y _{a}* and

*Y*, are associated with eigenvectors that oscillate at the same frequency but shifted

_{b}*π*/2 one eigenvector with respect to the other. Although they are not connected by their respective eigenvalues’ uncertainties, we could still generate a collection of frames by reconstructing a filtered version of the original data with only these two principal components. The resulting spatial-temporal structure is the “quasi-monochromatic process” of the kind proposed in the previous paragraph. This process usually represents a large amount of the total variance of the original data. The reconstruction of the original data with these two eigenvectors can be written as follows,

where the superindex ^{a,b} stands for the two principal components involved in the calculation. The eigenvectors, *e _{a}* and

*e*, are given by

_{b}Here the reconstructed frames can be given as

where we set *ϕ* = 0 without loss of generality. Equation (7) defines a complex principal component, *Y*̃_{a,b}, that deserves special attention. Alternatively, it can be written as follows,

where *α*(*x*,*y*) is the phase of the complex principal component. Then Eq. (7) becomes

The phase of this principal component shows temporal and spatial dependence. The previous development is based on the linear character of the Maxwell’s equations and the media involved in the calculation.

As illustrated above, the application of the PCA method to FDTD data, obtained for a well-referenced photonic crystal simulation that generates a monopolar mode, has produced a couple of eigenvectors with quasi-harmonic dependence. They can be grouped into a “quasi-monochromatic” process. From the results of the simulation this process represents 99.98701% of the total variance of the original data. In Fig. 4 (left), we show the results after reconstructing the frames with only the quasi-monochromatic process, *F*
^{1,2}(*x*,*y*, *t*), following Eq. (9). This movie can be seen as a filtered version of the original data obtained directly, without any other filtering technique, from the FDTD algorithm. The difference between the original and the PCA-filtered data is also included in this figure. It contains two main features. In the center, we can observe a high-frequency evolution. These patterns were also plotted in Fig. 2 and should be associated with noise processes. The other feature will be explained in the next subsection.

Figure 5 shows the individual eigenimages, *Y*
_{1}(*x*,*y*) and *Y*
_{2}(*x*,*y*). From these images, it is possible to obtain the modulus map and the phase map associated with *F*
^{1,2}(*x*,*y*, *t*). These maps may explain the existence and confinement of the electromagnetic radiation within the photonic crystal for the monopolar mode at its resonant frequency. The phase map has been obtained after applying an appropriate phase unwrapping algorithm [12]. The same procedure can be adopted for establishing the phase map of any mode of the microcavity, because the same structure of eigenimages appears systematically in any principal component decomposition of the electric field of any mode. When some other mode (hexapolar and quadrupolar modes) are excited in the photonic crystal, we have been able to define the corresponding “quasi-monochromatic”
processes that allow the calculation of the distribution associated with each one of the modes. The determination of the modulus and the phase related to these modes are presented in Fig. 6.

#### 3.2. Standing waves

Another interesting result is related to the characteristics of those principal components having quasi-harmonic dependence in their temporal evolution that can not be grouped with another principal component to form a “quasi-monochromatic” process. The associated eigenvectors are quasi-harmonic and are characterized by a given temporal frequency. No other principal component shows temporal evolution with quasi-harmonic dependence at the same frequency. Typically, the temporal frequency of the associated eigenvector falls in the tails of the spectrum for some of the modes of the photonic crystal. These structures are labelled as standing waves. In some other previous contribution it has been stated the presence of stationary waves within the structure of photonic crystals [13].

In the analyzed examples, the frequencies of the standing waves are located outside the band gap of the structure[6]. The band-gap comprises frequencies between $0.29\frac{c}{a}$ and $0.42\frac{c}{a}$, where *c* is the speed of light and *a* is the lattice constant of the crystal. We have simulated three modes: monopolar, one of the quadrupolars, and the hexapolar. The standing waves only appear distinguished from noise for the monopolar and the hexapolar modes computations. In the monopolar case, the frequency of the standing wave is $0.281\frac{c}{a}$, that is smaller than the lower limit of the bandgap. Moreover, for the hexapolar case, the frequency, $0.57\frac{c}{a}$, is higher than the upper limit of the bandgap. This is because the tails of the monopolar and hexapolar spectra spread outside the band gap with a non-negligible amplitude. This is not the case for the quadrupolar mode, whose tails outside the bandgap are negligible [6].

In order to assure that these processes are actually present in the spatial-temporal evolution of the data, we have performed several simulations by using excitation sources with the characteristic temporal frequency of the obtained eigenvector. The results have always produced only one principal component having such a temporal evolution. For the monopolar mode, the PCA has produced a principal component, *Y*
_{3}(*x*,*y*), that shows a harmonic dependence, *e*
_{3}(*t*) (see figure 3.b). The contribution of this third principal component to the total variance of the data is only *w*
_{3} = 0.00971%. Then, the presence of this oscillation may be swept out by any filtering technique applicable to the original data. On the other hand, the spatial distribution is quite remarkable and resembles the geometrical arrangements of the photonic crystal. Besides, the associated eigenvector follows a quasi-harmonic dependence. When reconstructing a filtered set of frames by using only this principal component, we may observe the typical behavior of an standing wave (see left map in Fig. 7). The nodes are also located in good accordance with the symmetry of the arrangement. The modulus of the instantaneous Poynting vector is displayed in the center of figure 7. This plot suggests that the energy is trapped at the cylinders’ locations. The hexapolar mode has produced another standing wave principal component, *Y*
_{3}(*x*,*y*) (see right map of figure 7). Its contribution to the total variance is 0.00569%. The symmetry of this standing wave is driven by the symmetry of the hexapolar mode and it is again attached to the
central cylinder structure.

#### 3.3. Decentering artifacts

The analysis of the computational noise generated by the FDTD algorithms is also possible. One of the most striking results of the usefulness of the PCA has been obtained when the monopolar mode of the analyzed photonic crystal is generated within a two dimensional grid with even number of nodes along the orthogonal directions. In such case, the excitation can not be precisely located at the center of the photonic crystal. The excitation frequency is that one of the monopolar mode. However, the application of the PCA has produced as a secondary process the excitation of noncentered modes, as the hexapolar one. The monopolar mode is still the most contributing one, but a detailed analysis of the obtained eigenimages has revealed the presence of another oscillation, clearly linked with the lack of centration of the excitation within the two-dimensional grid. This oscillation follows the spatial and temporal pattern of the hexapolar mode (see Fig. 8).

## 4. Conclusions

The use of FDTD techniques for simulating the interaction of electromagnetic fields with material structures has become a well-established tool for the analysis of photonic crystals. When PCA techniques are applied to the temporal evolution of the maps of electromagnetic magnitudes, the obtained results reveal spatial-temporal structures that can be interpreted and linked with the actual geometrical and constitutive properties of the photonic crystals. In this paper, we have demonstrated the pertinence of such analysis of FDTD results. The outputs of the PCA method have been interpreted as follows: the eigenimages represent maps of electromagnetic fields (here, the electric field and the Poynting vector); the eigenvectors represent temporal evolutions linked with the expected harmonic variations of the modes of the photonic crystal; and the eigenvalues represent their contribution to the total variance of the original data.

PCA allows an easy classification and grouping of the principal component into processes. When the classification of these processes is made in terms of eigenvalue distinguishability, the number of principal components associated with physically meaningful maps is drastically reduced, providing a simple tool for cleaning and filtering the FDTD results. Inasmuch as the PCA analyzes the correlation structure of the data, some spatial-temporal patterns with different correlations are distinguished. Furthermore, these structures with some correlation are extracted from the noise itself. When PCA is applied to the photonic crystal analysis, it has produced three types of spatial-temporal patterns. The first one is a “quasi-monochromatic” process related to the defect modes of the photonic crystal. It gives the resonant frequency, and the spatial distribution of its modulus and phase. The second one has been interpreted as a standing wave linked to the geometry and symmetry of the microcavity and the mode. Their frequencies fall outside the bandgap. Finally, the third spatial-temporal pattern can be interpreted as numerical noise. Some of the spatial-temporal structures with physical meaning contribute with much less than 1% to the total variance; however, their spatial shape and their temporal evolution denote inner relations that can be properly determined. In this case, they could be swept out by conventional filtering techniques.

When symmetric structures are simulated with symmetric calculation grids, special attention has to be paid to the matching of these symmetries. PCA has shown that even a minimum decentering may provoke the excitation of other modes of the photonic crystal. Owing to the negligible contribution to the total variance of that mode, it is not surprising that the application of any other filtering technique hides and erases the existence of those decentering artifacts. However, PCA has revealed them and quantified their importance. Moreover, some other subtle effects and interpretations are possible by a dedicated analysis of the data.

## Acknowledgments

This research has been supported by the project “Optical Antennas,” TIC2001-1259, financed by the Ministerio de Ciencia y Tecnología of Spain.

## References and links

**1. **E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. **58**, 2059–2062 (1987). [CrossRef] [PubMed]

**2. **P. R. Villeneuve, S. Fan, and J. D. Joannopoulos, “Microcavities in photonic crystals: mode symmetry, tunability, and coupling efficiency,” Phys. Rev. B **54**, 7837–7842 (1996). [CrossRef]

**3. **H. Y. Ryu and M. Notomi, “Enhancement of spontaneus emission from the resonant modes of photonic crystal slab single-defect cavity,” Opt. Lett. **28**, 2390–2392 (2003). [CrossRef] [PubMed]

**4. **C. Lixue, D. Xiaxou, D. Weiqiang, C. Liangcai, and L. Shutian, “Finite-difference time-domain analysis of optical bistability with low treshold in one-dimensional non-linear photonic crystal with Kerr medium,” Opt. Commun. **209**, 491–500 (2002). [CrossRef]

**5. **A. Taflove and S. C. Hagness*Computacional Electrodynamics: The Finite-Difference Time Domain Method*, 2nd ed. (Artech House, Boston, 2000).

**6. **S. Guo and S. Albin, “Numerical techniques for excitation and analysis of defect modes in photonic crystals,” Opt. Express **11**, 1080–1089 (2003); http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-9-1080 [CrossRef] [PubMed]

**7. **J. Jiang, J. Cai, G. P. Nordin, and L. Li, “Parallel microgenetic algorithm design for photonic crystal and waveguide structures,” Opt. Lett. **28**, 2381–2383 (2003). [CrossRef] [PubMed]

**8. **J. M. López-Alonso, J. Alda, and E. Bernabéu, “Principal component characterization of noise for infrared images,” Appl. Opt. **41**, 320–331 (2002). [CrossRef] [PubMed]

**9. **D. F. Morrison, *Multivariate Statistical Methods*, 3rd ed. (McGraw-Hill, Singapore, 1990), Chap. 8.

**10. **J. M. López-Alonso and J. Alda, “Operational parametrization of the 1/f noise of a sequence of frames by means of the principal components analysis in focal plane arrays,” Opt. Eng. **42**, 427–430 (2003). [CrossRef]

**11. **R. B. Cattell, “The scree test for the number of factors,” J. Multivar. Behav. Res. **1**, 245–276 (1966). [CrossRef]

**12. **D. D. Ghiglia and M. D. Pritt, *Two Dimensional Phase Unwrapping* (Wiley, New York, 1998), Chap. 4.

**13. **A. Tredicucci, “Marriage of two device concepts,” Science **302**, 1346–1347 (2003). [CrossRef] [PubMed]