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

The local density of states and response to an incident plane wave of a finite sized photonic crystal (PC) with nonlinear material (NLM) is analyzed. Of particular interest is the excitation of surface wave modes at the truncated surface of the PC, which is collocated with the NLM material. We compute the 2D Green function of the PC with linear material and then include the Kerr NLM in a self-consistent manner. The 2D PC consists of a square array of circular rods where one row of the rods is semi-circular in order to move the surface wave defect mode frequency into the band gap. Since the surface modes are resonant at the interface, the NLM should experience at least an order of magnitude increase in field intensity. This is a possible means of increasing the efficiency of the PC as a frequency conversion device.


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]. 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] In this paper, we compute LDOS and field characteristics for a finite two-dimensional PC consisting of both linear and nonlinear material. Some researchers [8] 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 Sec.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 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. Sec.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.

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.

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. [10,11,12] 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 . 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 [9]. dyadic G =ẑẑG zz effectively becomes a scalar that satisfies Dyson's equation and the electric field E =ẑE z satisfies either of the following forms of the Lippmann-Schwinger equation, where k 0 = ω/c and ρ = (x, y). The integration in Eqs. (1)- (3) is only over the cross-sectional area of the cylinders.
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 0 z (ρ) 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. [13] 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ε 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 To initiate the iteration algorithm, set n = 1 and there is a single perturbation ε k 1 at grid k 1 . The M 0 k 1 term represents the integration term when i = j = k 1 as where R = ∆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 where E 0 i is the incident field at the ith coordinate and N is the number of perturbation elementŝ ε k n in the structure. We have given a very brief outline of the method since details [10,11,12] are described elsewhere and the method is not the main focus of this work.

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 where G 0 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], the local density of states is written as LDOS = − 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 (ρ, ρ) = 0.25. It follows that the normalized LDOS given in numerical results below can easily be compared to free space relative to 0.25.

Non-linear material
Others have considered the effect of NLM on the infinite photonic crystal. [14,15,16,17,18] Using the iterative procedure [10,11,12] 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. [15] Assuming a Kerr non-linearity we writê ε (ρ) = ε c − 1 + χ|E c (ρ)| 2 : ρ within any semicircular cylinder ε c − 1 : ρ within any round cylinder 0 : ρ otherwise (9) Initially, the electric field distribution within the semi-circular rods, E 0 c (ρ), is obtained when ε 0 (ρ) = ε c − 1. Using E 0 c (ρ), the permittivity is then changed toε 1 (ρ) = ε c − 1 + χ|E 0 c (ρ)| 2 . Usingε 1 (ρ), a new Green function G i j and field distribution E 1 c (ρ) are calculated. This process 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| 3 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. 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 (n) (x, y) for all M points. We assume convergence is satisfied when 1 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.

Numerical Results
The PC is modeled as a square array of dielectric rods with period a and radius r with physical parameters [9] 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 ω 0 within the band gap where (ω 0 /c)(a/2π) = 0.350. We set N = 20 and this yields a free space resolution of ∆ = 0.0175λ 0 with ω 0 /c = 2π/λ 0 . The resolution in a cylinder is then degraded to ∆ c = ∆ √ 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 N zz (ρ, ρ ′ ). We find that the electric fields computed by both methods are in excellent agreement and we conclude that the method of calculating G N zz (ρ, ρ ′ ) is also accurate. 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 o , 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)| 3 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 o and 135 o , 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)| 3 results when χ = 0.0 and 0.01. The addition of χ = 0.01 causes a slight decrease in |E(x, y)| 3 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.

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.