Local density of states analysis of surface wave modes on truncated photonic crystal surfaces with nonlinear material

J. Merle Elson; Klaus Halterman

doi:10.1364/OPEX.12.004855

1. Introduction

The refinement of lithographic and other fabrication techniques used in photonic crystal (PC) development has generated a broad range of experimental, and theoretical interest [1, 2, 3, 4, 5, 6, 7, 8]. In particular, the ability of PC structures to control the propagation of electromagnetic fields has stimulated certain emerging concepts and devices that are based on PC technology, including microscopic lasers, optical switches, light localization, and resonant cavity antennas [4].

It is computationally convenient when investigating PC structures, to assume the PC is of infinite extent. This means a perfect crystal, or a crystal with defects in a supercell configuration, and this works well for band structure calculations. This approximation has limitations, since experimental situations involving finite-sized crystals may involve defect modes that result from truncating the infinite periodicity, or reflections and interference effects of the fields around the boundaries [2]. The vast majority of theoretical studies also assume the PC consists of linear material. When characterizing the electromagnetic properties of PC structures, the transmittance characteristics often yield limited information regarding the band gap, owing mainly to the directional dependence of the incident field [6]. The quantity of experimental relevance, that affords a more complete and fundamental picture of the spectra and modes of the system, is the spatially and energy resolved local density of states (LDOS) [7, 8].

In this paper, we compute LDOS and field characteristics for a finite two-dimensional PC consisting of both linear and nonlinear material. Some researchers [9, 10] have focused on surface modes associated with a one-dimensional PC. Our emphasis is on the properties of surface resonant modes as a localized field to concentrate the electric field in the presence of NLM. For the case of two-dimensional truncated structures consisting of an array of cylinders, surface electromagnetic waves can be supported. These surface modes are localized near the truncation boundary and propagate parallel to the interface. Since the intensity of the surface modes is concentrated at the interface, and if the interface contains nonlinear material (NLM), then a more efficient frequency conversion might be achieved. This could be useful in frequency conversion devices for use in detectors to operate between visible and infrared wavelengths.

In Section 2, we outline the integral formulation of the problem. Next, we briefly discuss the iterative numerical algorithm used to calculate the electromagnetic fields and Green’s function, which gives essential information regarding the LDOS. The LDOS will help reveal certain defect modes of the PC that occur at frequencies within a band gap. We present the results in Section 3 and find that within the PC structure, the LDOS is small but nonzero within the band gap, and that any localized states reside mainly within the truncated surface. Results are also shown that illustrate the effect of incorporating NLM material.

Fig. 1. Schematic of photonic crystal (PC) cross-section. The PC consists of infinitely long rods in a 2D array of 7 rods × 6 rods. Note that the rods in the topmost layer have been truncated along their diameter into a semi-circular shape. This truncation will allow localized surface modes to exist along this side of the PC within band gap frequencies. The lines AC and BD are guides for the eye in relating to numerical results below. The arrow and angle θ denote an incident plane wave. The PC has period a, the cylinders have radius r = 0.2a with permittivity ε_c = 8.9, and the background medium is vacuum. Dispersion characteristics of this PC are given in Fig. 2.

Download Full Size | PDF

2. Theory

We calculate the LDOS of a finite size 2D PC for the cases where the dielectric material is all linear or some rods are nonlinear. We will also investigate the response of the PC to an incident electric field. Of particular interest are localized surface wave modes along a PC interface. To model the 2D PC, we consider a finite-sized square array of dielectric cylinders as shown in Fig. 1. We begin with the case where the PC consists of linear material and assume the electric vector is polarized parallel to the rods. Since the PC considered here is finite size, the possibility of defect and localized modes exist. The design parameters of the PC are such that a band gap exists and truncation of the topmost layer of rods permits localized surface wave modes to occur at frequencies within the band gap. The dispersion of the bulk PC and surface modes is shown in Fig. 2.

In the case of linear media, the LDOS is independent of an incident field. In the case of non-linear material, we illuminate the PC with a plane wave and include non-linear material effects in a self-consistent manner. In this case, any change in the LDOS as calculated here, is examined. We neglect the effects of coupling to higher harmonic fields.

Fig. 2. (a) Band diagram of the bulk PC for an infinite 2D array of round rods with the electric vector parallel to the rods (E_∥) (blue) and the magnetic vector parallel to the rods (H_∥) (red). The vertical axis is normalized frequency f = (ω/c)(a/2π). As shown by the blue curves, an infinite PC of all round rods with the parameters given in Fig. 1 has a complete band gap in the frequency range 0.32 < f < 0.44 for (E_∥) modes. For (H_∥) modes, there is no complete bandgap. (b) These data are a subset of the Γ-X portion of (a) where the green curve shows the dispersion of E_∥ surface modes that propagate with wavevector k parallel to the semi-circular cylinder interface at frequencies within the bandgap. The black line is the light line and the gray lines show the band gap limits. Of particular interest in the analysis to follow is the green curve in (b) for the localized surface wave mode. We have chosen to truncate the top row of cylinders along their diameter, however varying the amount of truncation could be done and this would yield different dispersion curves for the surface mode. The data in (a) and (b) are reproduced from [11].

Download Full Size | PDF

2.1. Dyson’s equation and iterative solution

The method of analysis involves iterative solution to Dyson’s equation to calculate the Green function and the Lippman-Schwinger equation for the electric field [12, 13, 14]. We assume harmonic fields that vary as exp(-iωt). Here we consider an especially simple case where the electric vector is polarized parallel to the axis of the dielectric cylinders. In this case, the Green dyadic G = ẑẑG_zz effectively becomes a scalar that satisfies Dyson’s equation

G_{z z} (ρ, ρ′) = G_{z z}^{0} (ρ, ρ′) + k_{0}^{2} \int_{A} d^{2} ρ′′ G_{z z}^{0} (ρ, ρ′′) \hat{ε} (ρ′′) G_{z z} (ρ′′, ρ′)

and the electric field E = ẑE_z satisfies either of the following forms of the Lippmann-Schwinger equation,

E_{z} (ρ) = E_{z}^{0} (ρ) + k_{0}^{2} \int_{A} d^{2} ρ′ G_{z z} (ρ, ρ′) \hat{ε} (ρ′) E_{z}^{0} (ρ′)

E_{z} (ρ) = E_{z}^{0} (ρ) + k_{0}^{2} \int_{A} d^{2} ρ′ G_{z z}^{0} (ρ, ρ′) \hat{ε} (ρ′) E_{z} (ρ′)

where k ₀ = ω/c and ρ = (x,y). The integration in Eqs. (1)–(3) is only over the cross-sectional area of the cylinders. The $G_{z z}^{0} (ρ, ρ′) = \frac{i}{4} H_{0} (k_{0} ∣ ρ - ρ′ ∣)$ is the Green function for free space in the absence of the PC structure and G_zz is the Green function with the effects of the PC included. The $E_{z}^{0}$ (ρ) is the incident electric field. Note that when ρ = ρ́ in Eqs. (1)–(3) the singularity associated with the Green function would normally require special treatment, but in the simple case treated here, the singularity is more easily handled.[15] For simplicity, we omit all z subscripts henceforth.

The PC structure is described by the ε̂(ρ) that is non-zero only within the dielectric rods. Initially, we assume the rods are homogeneous with permittivity ε_c embedded in a background medium of unit permittivity. This yields

\hat{ε} (ρ) = ε_{c} - 1

when ρ is within any rod and zero otherwise. Numerical solution to Eqs. (1) and (3) can be performed by direct inversion, but this would typically be very large in terms of memory requirements, especially in Eq. (1). Here we provide only a brief description of the iterative method. The computational space is divided into an array of grids where the ith grid has area ΔA. The center of the ith and jth grid is located at ρ_i and ρ_j . The space is initially empty and the PC is then constructed by adding one “infinitesimal” ε̂ element at a time. When the nth element ε̂k_n is added we have

G_{k_{n} j}^{n} = {(I - M_{k_{n}}^{n - 1} {\hat{ε}}_{k_{n}})}^{- 1} G_{k_{n} j}^{n - 1} i = k_{n}, j \neq k_{n}

where $k_{0}^{2}$ ΔA = Â. To initiate the iteration algorithm, set n = 1 and there is a single perturbation ε _k₁ at grid k ₁. The M ⁰ _k₁ term represents the integration term when i = j = k ₁ as

M_{k_{1}}^{0} = k_{0}^{2} \int_{A} G_{k_{1} k_{1}}^{0} d A = \frac{i π k_{0} R H_{1} (k_{0} R)}{2} - 1

where $R = \sqrt{Δ A / π}$ . In (1), the integral over the grid square has been approximated by integrating over an equivalent circular area of radius R. The electric field solution can be built up in a similar iterative fashion or calculated as

E_{i} = E_{i}^{0} + \hat{A} \sum_{n = 1}^{N} G_{i k_{n}}^{N} {\hat{ε}}_{k_{n}} E_{k_{n}}^{0}

where $E_{i}^{0}$ is the incident field at the ith coordinate and N is the number of perturbation elements ε̂k_n in the structure. We have given a very brief outline of the method since details [12, 13, 14] are described elsewhere and the method is not the main focus of this work.

2.2. Local density of states

Computation of the electric field and the Green dyadic yields much information, including the normalized LDOS, that is defined here as

LDOS (ρ) = \frac{ℑ [G (ρ, ρ)]}{ℑ [G^{0} (ρ, ρ)]},

where G ⁰ is the free space Green function as given in Eq. (1), and following Eq. (3), and ℑ indicates the imaginary part. In the case of linear media, the LDOS as calculated from Eq. (8) is only dependent on the material and structural parameters of the PC. If nonlinear material is introduced, then obviously an excitation field is required to activate the nonlinearity. In this case, we still use Eq. (8) to calculate the LDOS with the PC illuminated by an incident plane wave of unit amplitude and at a given angle of incidence. This is done self-consistently yielding changes in the permittivity as described in the next section. The resulting changes in permittivity yields different values for G(ρ,ρ) and consequently the LDOS(ρ). The modeling of the nonlinearity is discussed in the following section.

In [7, 8], the local density of states is written as $LDOS = - \frac{2 ω}{π c^{2}} ℑ [G (ρ, ρ)]$ . Since in this work we consider the LDOS in Eq. (8) that is normalized to free space, we note that, -ℑG ⁰(ρ,ρ) = 0.25. It follows that the normalized LDOS given in numerical results below can easily be compared to free space relative to 0.25.

2.3. Non-linear material

Others have considered the effect of NLM on the infinite photonic crystal [16, 17, 18, 19, 20, 21]. Using the iterative procedure[12, 13, 14] briefly described above, we calculate for a finite size PC the Green function, LDOS, and the response of the PC shown in Fig. 1 to an incident plane wave. To consider the addition of a NLM in the topmost layer of semi-circular rods, we take a self-consistent approach [18]. Assuming a Kerr non-linearity we write

\hat{ε} (ρ) = {\begin{matrix} ε_{c} - 1 + χ {∣ E_{c} (ρ) ∣}^{2} : ρ within any semicircular cylinder \\ ε_{c} - 1 : ρ within any round cylinder \\ 0 : ρ otherwise \end{matrix}

Initially, the electric field distribution within the semi-circular rods, $E_{c}^{0}$ (ρ), is obtained when ε̂⁰(ρ) = ε_c - 1. Using $E_{c}^{0}$ (ρ), the permittivity is then changed to ε̂¹(ρ) = ε_c - 1 + χ| $E_{c}^{0}$ (ρ)|². Using ε̂¹(ρ), a new Green function G_ij and field distribution $E_{c}^{1}$ (ρ) are calculated. This process is repeated until succeeding electric fields throughout the PC converge to a stable result. This yields not only the electric fields, but also any changes in the LDOS.

In testing for convergence of the self-consistent algorithm, we found that for χ ≤ 0.01 approximately 10 iterations or less were needed to converge. If the discretized computational domain contains M points (x,y), we calculate the electric field E ⁽ⁿ⁾(x,y) for all M points. We assume convergence is satisfied when $\frac{1}{M} \sum_{x, y} ∣ r^{(n)} (x, y) ∣ \leq 0.0001$ where r ⁽ⁿ⁾(x,y) = [E ⁽ⁿ⁾(x,y) - E ^(n-1)(x,y)]/E ^(n-1)(x,y) and n is the number of iterations. We did not do an extensive parametric evaluation of the self-consistent convergence conditions. We did note that when χ = 0.1, the self-consistent algorithm failed to converge. In any case, we did not try to use experimentally obtained χ values, rather we note that χ ≤ 0.01 is well within any realistic material values.

3. Numerical results

The PC is modeled as a square array of dielectric rods with period a and radius r with physical parameters [11] chosen so that there is a complete band gap for the infinite crystal and this is shown in Fig. 2(a) and the defect mode surface wave dispersion is shown in Fig. 2(b). The computational space is divided into square grids of side dimension Δ and the PC period is a = NΔ. We choose a frequency ω ₀ within the band gap where (ω ₀/c)(a/2π) = 0.350. We set N = 20 and this yields a free space resolution of Δ = 0.0175λ ₀ with ω ₀/c = 2π/λ ₀. The resolution in a cylinder is then degraded to $Δ_{c} = Δ \sqrt{8.9}$ , and this is still adequate for accurate computation.

To check the accuracy of our results, we have compared calculation of the electric field as obtained from Eqs. (2) and (3). We note that for reasonable sized grid of N points, Eq. (2) can be solved for E_z (ρ) by direct inversion. Also, E_z (ρ) can be obtained from Eq. (3) after first obtaining $G_{zz}^{N}$ (ρ,ρ́). We find that the electric fields computed by both methods are in excellent agreement and we conclude that the method of calculating $G_{zz}^{N}$ (ρ,ρ́) is also accurate.

Fig. 3. (a) Normalized LDOS(x) vs position x within the PC for various frequencies with linear material and (b) Normalized LDOS(x) at resonance frequency compared with nonlinear material. In (a), and referring to Fig. 1, the LDOS(x) is shown on a linear scale and is calculated along line AC (between rods) and BD (through rods) where the BD data has been displaced one unit upward for clarity. For reference, the large black circles in the lower part of the figures indicate the positions of the rods in relation to the LDOS. In (a), the dimensionless frequency parameter f = (ω/c)(a/2π) is varied for ω within the band gap from 0.34 to 0.36 with a peak at f = 0.344. Along line BD, there is clear indication of localized states concentrated at the semi-circular rod interface (x ≈ 130Δ, where Δ is the spatial resolution) whereas no localization at the opposite circular rod interface (x ≈ 10Δ). Along line AC there is no comparable indication of localized states. It follows that the localized states exclusively reside within the semi-circular rods. This indicates the importance of defect parameters, such as altering the rod cross-section, in the surface mode dispersion. In (b), the resonant LDOS(x) is plotted on a log scale along lines AC and BD for f = 0.344 with χ = 0 and 0.01. The χ = 0 curves are repeated from (a). To activate the NLM, the PC is illuminated by a plane wave of unit amplitude incident at θ = 45°. Along line AC, the green and black curves are superimposed. Along line BD, the red and blue curves are nearly superimposed with the blue curve slightly larger in value that indicates a slight reduction in LDOS with non-zero χ. This reduction is not very clear on this log scale, however the log scale does show more clearly the LDOS(x) throughout the PC. This indicates that the presence of non-linear material in the semi-circular rods does not significantly alter the LDOS. Finally, these numerical data are completely independent of an incident field.

Download Full Size | PDF

Fig. 4. Plane wave at normalized frequency f = 0.344 incident on PC at angles of (a) 45° and (b) 135° (see Fig. 1) for both linear and non-linear material. The x-coordinate is in Δ units and the electric field is plotted on a log scale as |E|³ since this is relevant to a Kerr nonlinearity. For both angles of incidence, there is significant excitation of the surface defect mode as noted by the peaks around x ≈ 130Δ and lines BD. In (a), there are incident and reflected fields in the x ≤ 0 PC interface region. In (b), there are again the incident and reflected fields for x ≥ 130Δ in addition to the surface mode excitation field. Note also in (b) that the low field intensity for x ≈ 0 region is because this is a near-field shadow region. It is seen that for a unit amplitude incident wave, the localized |E|³ field is well in excess of an order of magnitude larger because of the surface mode excitation and the localization is precisely at the sites of the nonlinear material.

Download Full Size | PDF

Fig. 5. Contour plots of a plane wave at frequency f = 0.344 incident on PC at angles of (a) 45° and (b) 135° for linear material. The angles of incidence are indicated by the arrows. These plots show the electric field intensity throughout the PC and within 2 periods a outside the PC., The line plots of Fig. 4 are a subset of these data. The x and y-coordinates are in λ ₀ = a/0.35 units. Surface defect mode excitation is clearly seen along the truncated layer of rods. Note that the outline of the PC rod structure is superimposed on this contour plot.

Download Full Size | PDF

Figure 3 shows the results of computations of the LDOS associated with the PC and the plane wave response of the PC. In Fig. 3(a), the calculated LDOS is shown for several frequencies and a PC with all linear material. At the surface mode resonant frequency, and with the PC illuminated with a unit-amplitude plane wave at θ = 45°, the computed effect of including NLM is shown in Fig. 3(b). The NLM coefficient is chosen as χ = 0.01. The results of the LDOS calculations indicate a large concentration of states available at the truncated-rod interface. When the NLM coefficient χ is non-zero, there is noticeable decrease in the LDOS. This decrease is due to the change in the permittivity of the semi-circular rods and the surface wave mode is quite sensitive to any physical or material changes in the semi-circular interface.

In Fig. 4 we consider the response of the PC when illuminated by a plane wave at two different angles relative to the truncated rod interface. We plot the electric field as |E(x,y)|³ since this is the magnitude relevant to Kerr nonlinearity. For both angles of incidence, there is strong excitation of the surface mode. It is seen that for χ = 0.0 and 0.01 for both angles of incidence θ = 45° and 135°, the surface mode peaks are quite similar in both cases. This is likely due to the fact that very little of the PC actually contains nonlinear material. Nevertheless, there are some differences in the |E(x,y)|³ results when χ = 0.0 and 0.01. The addition of χ = 0.01 causes a slight decrease in |E(x,y)|³ and this is consistent with the LDOS data in Fig. 4(b).

A more comprehensive view of the field profile is illustrated in Fig. 5, which shows the corresponding contour plots relating to Fig. 4 for linear media. Here the field magnitude is plotted over the entire PC including two periods outside the PC boundary. The surface mode excitation at the semicircular rod sites is clearly seen. The field magnitude inside the PC is quite diminished and the location of the rods comprising the PC is shown by an overlay of an outline of the rods.

4. Conclusions

An analysis of the LDOS for surface wave modes on a truncated 2D photonic crystal interface has been accomplished using a Green function formalism. The PC is composed of circular rods (except for the semicircular interface) and the electric field is parallel to the rods. The results show that the presence of NLM does not significantly alter the LDOS as defined in this paper and this could be due to the NLM being a small fraction of the total volume of PC material. If the idea of using a 2D PC as a component in frequency conversion is viable, then it seems that additional volume of NLM would probably be needed. We have also neglected the fact that frequency conversion does occur and that there is an additional coupled equation to consider that would describe the generation of higher harmonic fields due to the Kerr nonlinearity. This may be a source of damping that could broaden the resonance associated with the surface wave mode. This could reduce the conversion efficiency.

Acknowledgments

Support was provided by ONR Independent Laboratory In-house Research funds and DARPA MORPH funds. Valuable computing resources and assistance were available from the Arctic Region Supercomputing Center.

References and links

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

2. W. M. Robertson, G. Arjavalingam, R.D. Meade, K.D. Brommer, A.M. Rappe, and J.D. Joannopoulos, “Observation of surface photons on periodic dielectric arrays,” Opt. Lett. 18, 528 (1993). [CrossRef] [PubMed]

3. R. Hillebrand, St. Senz, W. Hergert, and U. Gösele, “Macroporous-silicon-based three-dimensional photonic crystal with a large complete band gap,” J. Appl. Phys. 94, 2758 (2003). [CrossRef]

4. 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 (2000). [CrossRef]

5. I. V. Konoplev, A. D. R. Phelps, A. W. Cross, and K. Ronald, “Experimental studies of the influence of distributed power losses on the transparency of two-dimensional surface photonic band-gap structures,” Phys. Rev. E 68, 066613 (2003). [CrossRef]

6. Y. Wang, B. Cheng, and D. Zhang, “Distribution of density of photonic states in amorphous photonic materials,” Appl. Phys. Lett. 83, 2100 (2003) ; Y.W. Wang, X. Hu, X. Xu, B. Cheng, and D. Zhang, “Localized modes in defect-free dodecagonal quasiperiodic photonic crystals,” Phys. Rev. B 68, 165106 (2003). [CrossRef]

7. D. P. Fussell, R. C. McPhedran, C. M. de Sterke, and A. A. Asatryan,“Three-dimensional local density of states in a finite-sized two-dimensional photonic crystal composed of cylinders,” Phys. Rev. E 67, 045601 (2003) [CrossRef]

8. A. A. Asatryan, K. Busch, R. C. McPhedran, L. C. Botten, C. M. de Sterke, and N. A. Nicorovici, “Two-dimensional Green’s function and local density of states in photonic crystals consisting of a finite number of cylinders of infinite length,” Phys. Rev. E 63, 046612 (2001). [CrossRef]

9. F. Villa, L. Regalado, F. Ramos-Mendieta, J. Gaspar-Armenta, and T. Lopez-Rios, “Photonic crystal sensor based on surface waves for thin-film characterization,” Opt. Lett. 27646 (2002) [CrossRef]

10. J. Gaspar-Armenta and F. Villa, “Photonic surface-wave excitation: photonic crystal-metal interface,” J. Opt. Soc. Am. B 202349 (2003) [CrossRef]

11. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, Princeton NJ, 1995), See page 75 for a diagram of the surface mode dispersion.

12. O. J. Martin, A. Dereux, and C. Girard, “Generalized field propagator for arbitrary finites-size photonic band gap structures,” J. Opt. Soc. Am. A 111073 (1994) [CrossRef]

13. O. J. Martin, C. Girard, D. R. Smith, and S. Schultz, “Iterative scheme for computing exactly the total field propagating in dielectric structures of arbitrary shape,” Phys. Rev. Lett. 82315 (1999) [CrossRef]

14. O. J. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E 58, 3903 (1998) [CrossRef]

15. A. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE 68, 248 (1980). [CrossRef]

16. P. Tran, “Photonic-band-structure calculation of material possessing Kerr nonlinearity,” Phys. Rev. B 52, 10673 (1995) ; [CrossRef]

17. P. Tran, “Photonic band structure calculation of system possessing Kerr nonlinearity,” in “Photonic Band Gap Materials” NATO ASI Series E: Applied Sciences 315, 555, Ed. C. M. Soukoulis (Kluwer Academic Publishers, Boston, 1996).

18. V. Lousse and P. Vigneron, “Self-consistent photonic band structure of dielectric superlattices containing nonlinear optical materials,” Phys. Rev. E 63027602 (2001) [CrossRef]

19. M. Soljačič, D. Luo, J.D. Joannopoulos, and S. Fan “Non-linear photonic crystal microdevices for optical integration,” Opt. Lett. 28637 (2003) [CrossRef]

20. M. Bahl, N. Panoiu, and R. Osgood, “Nonlinear optical effects in a two-dimensional photonic crystal containing one-dimensional Kerr defects,” Phys. Rev. E 67056604 (2003) [CrossRef]

21. S. Mingaleev and Y. Kivshar, “Nonlinear transmission and light localization in photonic-crystal waveguides,” J. Opt. Soc. Am. B 192241 (2002) [CrossRef]

Local density of states analysis of surface wave modes on truncated photonic crystal surfaces with nonlinear material

Abstract

1. Introduction

2. Theory

2.1. Dyson’s equation and iterative solution

2.2. Local density of states

2.3. Non-linear material

3. Numerical results

4. Conclusions

Acknowledgments

References and links

Cited By

Figures (5)

Equations (11)

Optics Express