1. Introduction

Photonic Crystals (PC’s) are periodic arrangements of dielectric or metallic materials. In these last years, dielectric PC’s have been intensively studied to realize optical devices, such as filters, resonators and waveguides, exploiting their ability to control light propagation [1–2]. It is well known that the unperturbed PC’s can exhibit forbidden bands, i.e. frequency ranges where no electromagnetic (e.m.) field is allowed. A single lattice defect can support a localized e.m. field within the frequency band gap of the unperturbed PC and a line of single defects can guide an optical signal along its direction creating frequency-selective waveguides. The importance of fast and accurate software tools for PC’s analysis and design is then evident. Among the different techniques, the plane-wave method (PWM) [1–5] and the Finite Difference Time Domain (FDTD) approach [6–8] are probably the most used ones.

The PWM, based on the Fourier expansion of the periodic dielectric constant, allows straightforward analysis of perfect PC’s. When defects are present, the PW approach can still be applied, but after the introduction of the so called supercell [1–4], a domain containing the defect and a part of the regular structure large enough to guarantee that the field vanishes at its boundary. It is then possible again to apply the PW expansion to a new structure made by periodically arranged supercells. In this case the complexity of the algorithm increases, in particular in the expression of the dielectric function, even if one can take advantage of the Fast Fourier Transform (FFT) [2] and of the shift property of the Fourier transform [5] to reduce CPU times. Particular care must also be used to choose the size of the supercell and the finite number of plane waves approximating the infinite sum as a compromise guaranteeing both convergence to correct results (large size and number of terms desired) and CPU time reduction (small size and few terms desired).

In this paper we present (Section 2) a new approach to the PWM, which reformulates the expression of the dielectric constant through the superposition of two regular lattice, the former for the perfect structure and the latter for the defects. This method leads to a single expression of the dielectric constant, which can be inserted in the eigenvalue and eigenmodes calculation routine. In Section 3, the method will be used to study a waveguide (W1) formed by single line defects in a two-dimensional triangular lattice of air holes in a dielectric substrate and will be shown to be able to accurately determine dispersion curves and field patterns in presence of “mini-stop bands”, evidencing also coupling phenomena [9–13] among fundamental and higher order modes (the so called curve “anticrossing”). These results will be compared also with those obtained by a 2D-FDTD method, to prove the correctness of our approach. Finally (Section 4), conclusions will be drawn.

2. Reformulation of the plane wave method

The method will be described, for simplicity, considering a 2D PC with square or triangular lattice, but it can easily be extended to any lattice shape and to 3D structures. We consider first the ideal 2D PC made by cylinders with square or triangular lattice, that we call the “perfect lattice”. The rods have their axis parallel to the z direction, lattice constant a, radius R and dielectric constant ε_a. The dielectric constant of the background material is ε_b. The PC is characterized by the lattice vectors {a₁,a₂} and the reciprocal lattice vectors {b₁,b₂} [1,2].

The relative dielectric constant of the structure without defects (denoted by subscript g), expanded in Fourier series in the reciprocal lattice, is

being r is the spatial vector position and G=h₁ b₁+h₂ b₂, with h₁,h₂ arbitrary integer numbers. The coefficient ε_g(G) is defined as

for G=0, while, for G≠0, it holds:

where $f_g=\frac{\pi R^2}{\mid {\mathit{a}}_{\mathit{1}}\times {\mathit{a}}_{\mathit{2}}\mid }$ is the filling factor of the perfect structure and J ₁ is the first order Bessel function.

Once the perfect lattice has been defined, we consider the presence of defects, corresponding, for example, to removal of rods. We then consider a new lattice formed by the defects repeated with proper periodicity along the two lattice directions, if single defects are present, or along one, if straight waveguides are considered. We call this new structure “defect perfect lattice”. The defect dielectric constant is ε_d=ε_b-ε_a while the background material has a dielectric constant equal to 0. The defect perfect lattice can be characterized by the vectors c₁=n₁ a₁ and c₂=n₂ a₂, with n₁, n₂ arbitrary integer numbers. The reciprocal lattice vectors are d₁=b₁/n₁, d₂=b₂/n₂. The generic vector in the reciprocal lattice is now defined as S=l₁d₁+l₂d₂, where l₁,l₂ are arbitrary integer numbers. The relative dielectric constant of the defect perfect lattice, denoted by subscript s, can then be expressed as

where the coefficient ε_s (S) is then

for S=0, while, for S≠0, it holds

being $f_s=\frac{\pi R^{\mathit{2}}}{\mid {\mathit{c}}_{\mathit{1}}\times {\mathit{c}}_{\mathit{2}}\mid }=\frac{f_g}{n_{\mathit{1}}{n}_{\mathit{2}}}$ the filling factor of the defects.

The PC with defects is then formed by the superimposed perfect lattice and defect perfect lattice. For example, considering a W1 waveguide, Fig. 1(a) shows the direct lattice, with defects painted in green, built with n₁=1 and n₂=3. Figure 1(b) reports the superposition of the two reciprocal lattices, the former (larger gray circles) referring to the perfect lattice, the latter (smaller green dots) referring to the defect perfect lattice. Note that the nodes of the perfect reciprocal lattice are a subset of the defect perfect reciprocal lattice.

The dielectric function formed by the perfect lattice and the defect perfect lattice, results adding the dielectric constant distributions (1) and (4) obtaining:

where ε_c(S) is the coefficient that describes the superposition of two reciprocal lattice. For nodes belonging to both reciprocal lattices, ε_c(S) can be written

for S=0, while, for S≠0

where l₁=m n₁ and l₂=n n₂, with m, n arbitrary integers different from 0.

 

Fig. 1. Modeling of 2D triangular lattice PC waveguide W1: (a) direct lattice of the W1 and (b) reciprocal lattices of the perfect lattice (larger gray circles) and the defect perfect lattice (smaller green dots).

Download Full Size | PDF

Otherwise, for nodes corresponding only to defects:

The structure usually obtained repeating the supercell is now built superposing two infinite periodic structures. This allows to calculate the dielectric constant distribution of the studied structure with a much simpler routine. Care must always be taken to choose the proper distance between the defects to avoid coupling phenomena, but this simply reduces to the best choice of n₁ and n₂, without the need to repeat the supercell construction procedure.

For the structure formed by the perfect lattice and the defect perfect lattice, the electromagnetic problem to be solved to determine eigenvalues and eigenmodes is the classical one [1,2] with the only formal difference that S is used instead of G:

where ω is the eigen-angular frequency and k the wave vector while the eigenfunctions of the electric and magnetic fields can be expanded in Fourier series [1,2]:

Note that the differences between this approach and the classical one, based on the supercell concept, refer only to the way to obtain the final Eqs. (11–12), which is now simpler and faster to run, and not to the results, which are always the same. We have tested the algorithm for perfect structures, single defect resonant cavities and waveguides and it proved stable and rapidly converging.

3. Results and comparison between different methods

As an example of application of the method, we consider a two-dimensional single line photonic crystal waveguide (W1) of air circular holes in a dielectric substrate, with relative dielectric constant equal to 11.56. The crystal has a triangular lattice and a filling factor of 83%. This structure has a photonic band gap both for TE and TM polarized fields. In the following we will restrict our attention to TM polarized fields (E parallel to the rod axes) propagating along the k direction (in the x-y plane). The method is applied considering waveguides repeated every 6 lattice periods a and 1089 plane waves. In this case, few minutes of CPU time are enough to run the code on a 1800 MHz Pentium IV Computer.

 

Fig. 2. (a) Dispersion curves computed by the PW method and (b) transmission coefficient computed by FDTD of a W1 waveguide. (c) Local picture of the “mini-stop band” associated to the coupling among the fundamental mode and higher order modes. (d)-(e)-(f) Electric displacement pattern associated to the fundamental mode. (g)-(h) Electric displacement pattern associated to the high order modes.

Download Full Size | PDF

The waveguide dispersion curves obtained with our method are shown in Fig. 2(a). The dispersion curves are calculated in the first Brillouin zone, in the y direction of Fig. 1, varying k between 0 and π/a. The gray areas are the projected band structure of the perfect crystal that limit the band gap between 0.46 a/λ and 0.54 a/λ. The solid lines are the dispersion curves of the waveguide modes. The red line, marked by index 1, refers to the fundamental mode.

Figure 2(a) shows that the dispersion curve of the fundamental mode crosses the dispersion curves of other modes. In particular, around 0.425 a/λ and 0.49 a/λ, the so called “anticrossing” phenomenon occurs, involving higher order modes. This generates two “mini-stop bands” [9–13], evidenced by the two circles, always in Fig. 2(a). Considering, in particular, the “mini-stop band” around 0.49 a/λ, Fig. 2(c) shows a close-up of the mini gap caused by the coupling between the fundamental mode and two higher order modes which are represented by blue and green lines, with index 2 and 3, respectively. Mode 3 is odd, while modes 1 and 2 are even. The field patterns of the waveguide modes corresponding to the five different points, marked in Fig. 2(c), are shown in the Figs. 2(d), 2(e), 2(f), 2(g), 2(h). In particular, Figs. 2(d), 2(e), 2(f) refer to the electric displacement of the fundamental mode. In the intermediate region (field of Fig. 2(e)), coupling of modes 1, 2 and 3 changes the field pattern, which results in a kind of zig-zag path caused by the interference of the odd and even modes. The even pattern of mode 2 is shown in Fig. 2(g), while the odd pattern of mode 3 is shown in Fig. 2(h).

 

Fig. 3. (1.7 MB) Movies of Electric field at (a) 0.47 a/λ and (b) 0.49 a/λ. [Media 2]

Download Full Size | PDF

The presence of these two “mini-stop bands” has been confirmed using a 2D-FDTD code, based on the Yee’s algorithm with total field-scattered field approach [6–8]. Figure 2(b) reports the calculated transmission characteristic of the structure. The whole FDTD computational domain is divided into 372×623 cells with step of Δx=Δy=0.005µm and 8 GT-PML layers. Agreement between the two set of results is very good. Finally Fig. 3(a) and 3(b) show how the 2D-FDTD accounts for propagation of the electric field in the structure at 0.47 a/λ and 0.49 a/λ, which are respectively in a allowed and in a “mini-stop band”: the attenuation of the field in the latter case is evident.

4. Conclusions

In this paper we have presented a new approach of the plane-wave method for photonic crystal analysis. This procedure expresses the dielectric constant distribution as the superposition of two functions, one for the perfect lattice and one for the defect perfect lattice. This leads to a simpler expression of the dielectric constant to be inserted in the standard eigenvalue problem, thus reducing the code overall complexity. The approach is easy to implement and fast to run. As an example, a 2D triangular lattice photonic crystal waveguide was studied, determining its dispersion curves and the guided field distributions, allowing also a detailed study of mode coupling effects and “mini-stop bands” presence. These results have been successfully validated by comparison with those obtained using a FDTD approach.