Coordinate transformation formulation of electromagnetic scattering from imperfectly periodic surfaces

Koki Watanabe; Jaromír Pištora; Yoshimasa Nakatake

doi:10.1364/OE.20.009978

1. Introduction

The periodic structures are widely used for optical, photonics, and electronics applications. The most representative example of the periodic structures may be the diffraction grating ruled on a plane substrate, and it plays important roles in spectroscopic measurements or optical communications. Recent improvement of the technology requires highly accurate numerical modeling in some cases. However, the real structures of gratings are not perfectly periodic due to the finite extent and fabrication error, and specifying the influence of the imperfection is sometimes important [1]. Besides, small collapses of the periodicity are sometimes useful to control the electromagnetic wave propagation in the periodic structures. Defects in the photonic crystals work as resonators or waveguides [2], and modulations of the local period in the gratings yield chirped gratings or the Fresnel lenses [3].

The aim of this paper is to present a rigorous spectral-domain approach to the electromagnetic scattering from periodically corrugated surface in which the structural periodicity is locally collapsed. The present formulation is based on the differential method of Chandezon et al. (C-method) [4], which is a well-known approach to the plane-wave scattering from surface-relief grating. When the plane-wave illuminates a perfectly periodic structure, the Floquet theorem asserts that the scattered fields are pseudo-periodic (namely, each field component is given by a product of a periodic function and an exponential phase factor). This implies that the fields have equidistant discrete spectra in the wavenumber space and can be expressed in the generalized Fourier series expansions, and many numerical approaches to periodic structures, including the C-method, utilize this property [5]. However, the structure considered in this paper is not perfectly periodic and the fields generally have continuous spectra. The spectral-domain approach introduces a discretization in the wavenumber space, and the present approach uses the pseudo-periodic Fourier transform (PPFT) [6] to consider the discretization scheme. The PPFT makes it possible to consider the scheme only inside the Brillouin zone because the transformed fields have a periodic property in terms of the transform parameter related to the wavenumber. Also, the transformed fields are pseudo-periodic in terms of the spatial parameters, and the conventional approaches to the periodic structures become possible to be applied for the scattering problem of imperfectly periodic structures. In Ref. [7], the scattering problem of the lamellar grating with defects is formulated by the rigorous coupled-wave analysis (RCWA) [8, 9] with the help of PPFT. RCWA is certainly suitable to analyze the lamellar gratings with deep depth. However, it is not always easy to apply for arbitrary shaped gratings [10, 11]. This paper combines the PPFT with the C-method and represents a part of wide applicability.

2. Settings of the problem

The present paper considers the scattering problem of electromagnetic fields with a time-dependence exp(−iωt) from imperfectly periodic surfaces, in which structural periodicity is locally collapsed. Figure 1 shows an example of the structures under consideration. The structure is uniform in the z-direction and the y-axis is perpendicular to the periodicity direction though the periodicity is not perfect. The equation of the corrugated surface is given by y = g(x), where g(x) is supposed to be a continuous function. The minimum and the maximum values of g(x) are denoted by h₁ and h₂, respectively. The surface profile function g(x) is supposed to be decomposed into periodic and aperiodic parts:

g (x) = g^{(p)} (x) + g^{(a)} (x)

where g^(p)(x) is a periodic function with the period d and g^(a)(x) has nonzero value only at a₁ ≤ x ≤ a₂ where the structural periodicity is locally collapsed. The surrounding region y > g(x) is filled with a homogeneous and isotropic medium with the permittivity ε_s and the permeability μ_s, and the substrate region y < g(x) is also filled with a homogeneous and isotropic medium described by the permittivity ε_c and the permeability μ_c. For lossy media, we use complex values for the permittivity and/or the permeability, and the terms concerning to the current densities are eliminated in the following formulation. The regions y > g(x) and y < g(x) are specified by s and c, respectively, and the wavenumber in each region is denoted by k_r = ω(ε_r μ_r)^1/2 for r = s, c. The electromagnetic fields are uniform in the z-direction and two-dimensional scattering problem is considered here. Two fundamental polarizations are expressed by TE and TM, in which the electric and the magnetic fields are respectively parallel to the z-axis. We sometimes denote the z-component of electric field for the TE-polarization and the z-component of magnetic field for the TM-polarization by ψ(x,y) to express both polarizations simultaneously. The incident field is supposed to illuminate the corrugated surface from the upper or lower regions and there exists no source inside the inhomogeneous region h₁ ≤ y ≤ h₂.

Fig. 1 An example of imperfectly periodic surfaces (surface-relief grating with a defect).

Download Full Size | PDF

3. Formulation

First, we consider the TE-polarized fields. From Maxwell’s curl equations and the constitutive relations, the Cartesian components of electromagnetic fields are described by the following relations:

\frac{\partial}{\partial x} H_{y} (x, y) - \frac{\partial}{\partial y} H_{x} (x, y) = - i ω ε_{r} E_{z} (x, y)

\frac{\partial}{\partial y} E_{z} (x, y) = i ω μ_{r} H_{x} (x, y)

\frac{\partial}{\partial x} E_{z} (x, y) = - i ω μ_{r} H_{y} (x, y)

in the region r (r = s,c). We introduce a curvilinear coordinate system O-uvz, which is related to the original coordinate system O-xyz by continuous transformation equations: u = x and v = y – g(x). Then, Eqs. (2)–(4) are transformed as follows:

[\frac{\partial}{\partial u} - ({\dot{g}}^{(p)} (u) + {\dot{g}}^{(a)} (u)) \frac{\partial}{\partial v}] H_{y} (u, v) - \frac{\partial}{\partial v} H_{x} (u, v) = - i ω ε_{r} E_{z} (u, v)

\frac{\partial}{\partial v} E_{z} (u, v) = i ω μ_{r} H_{x} (u, v)

[\frac{\partial}{\partial u} - ({\dot{g}}^{(p)} (u) + {\dot{g}}^{(a)} (u)) \frac{\partial}{\partial v}] E_{z} (u, v) = - i ω μ_{r} H_{y} (u, v)

where ġ^(p)(u) and ġ^(a)(u) denote the derivatives of g^(p)(u) and g^(a)(u), respectively.

Here, we introduce the PPFT. Let f(u) be a function of u and d be a positive real constant. Then the PPFT and its inverse transform are formally defined by

\bar{f} (u; ξ) = \sum_{m = - \infty}^{\infty} f (u - m d) e^{imd ξ}

f (u) = \frac{1}{k_{d}} \int_{- k_{d} / 2}^{k_{d} / 2} \bar{f} (u; ξ) d ξ

where ξ is the transform parameter and k_d = 2π/d. The transformed functions have pseudo-periodic property in terms of u: f̄(u − md;ξ) = f̄(u;ξ) exp(−imdξ) for any integer m, and also have periodic property in terms of ξ: f̄(u;ξ − mk_d) = f̄(u;ξ). We use the period of g^(p)(x) for the positive real constant d for the PPFT. Then, k_d becomes the inverse lattice constant and the periodicity cell of the transformed function gives the Brillouin zone. Applying the PPFT to Eqs. (5)–(7), we obtain the following relations:

\begin{array}{r} (\frac{\partial}{\partial u} - {\dot{g}}^{(p)} (u) \frac{\partial}{\partial v}) {\bar{H}}_{y} (u; ξ, v) - \frac{1}{k_{d}} \int_{- k_{d} / 2}^{k_{d} / 2} {\bar{\dot{g}}}^{(a)} (u; ξ - ξ^{'}) \frac{\partial}{\partial v} {\bar{H}}_{y} (u; ξ^{'}, v) d ξ^{'} \\ - \frac{\partial}{\partial v} {\bar{H}}_{x} (u; ξ, v) = - i ω ε_{r} {\bar{E}}_{z} (u; ξ, v) \end{array}

\frac{\partial}{\partial v} {\bar{E}}_{z} (u; ξ, v) = i ω μ_{r} {\bar{H}}_{x} (u; ξ, v)

\begin{array}{r} (\frac{\partial}{\partial u} - {\dot{g}}^{(p)} (u) \frac{\partial}{\partial v}) {\bar{E}}_{z} (u; ξ, v) - \frac{1}{k_{d}} \int_{- k_{d} / 2}^{k_{d} / 2} {\bar{\dot{g}}}^{(a)} (u; ξ - ξ^{'}) \frac{\partial}{\partial v} {\bar{E}}_{z} (u; ξ^{'}, v) d ξ^{'} \\ = - i ω μ_{r} {\bar{H}}_{y} (u; ξ, v) . \end{array}

Since the transformed fields are pseudo-periodic in terms of u, they can be approximately expanded in the truncated generalized Fourier series. For example, the z-component of electric field is written as

{\bar{E}}_{z} (u; ξ, v) = \sum_{n = - N}^{N} {\bar{E}}_{z, n} (ξ, v) e^{i α_{n} (ξ) u}

with

α_{n} (ξ) = ξ + n k_{d}

where N denotes the truncation order and Ē_z,n(ξ,v) are the nth-order coefficients. To treat the coefficients systematically, we introduce (2N + 1) × 1 column matrices; for example, the coefficients of Ē_z(u;ξ,v) are expressed by a column matrix ē_z(ξ,v) in such a way that its nth-entries are given by Ē_z,n(ξ,v). Then Eqs. (10)–(12) yield the following relations:

\begin{array}{r} (i \bar{U} (ξ) - 〚 {\dot{g}}^{(p)} 〛 \frac{\partial}{\partial v}) {\bar{h}}_{y} (ξ, v) - \frac{1}{k_{d}} \int_{- k_{d} / 2}^{k_{d} / 2} 〚 {\bar{\dot{g}}}^{(a)} 〛 (ξ - ξ^{'}) \frac{\partial}{\partial v} {\bar{h}}_{y} (ξ^{'}, v) d ξ^{'} \\ - \frac{\partial}{\partial v} {\bar{h}}_{x} (ξ, v) = - i ω ε_{r} {\bar{e}}_{z} (ξ, v) \end{array}

\frac{\partial}{\partial v} {\bar{e}}_{z} (ξ, v) = i ω μ_{r} {\bar{h}}_{x} (ξ, v)

\begin{array}{r} (i \bar{U} (ξ) - 〚 {\dot{g}}^{(p)} 〛 \frac{\partial}{\partial v}) {\bar{e}}_{z} (ξ, v) - \frac{1}{k_{d}} \int_{- k_{d} / 2}^{k_{d} / 2} 〚 {\bar{\dot{g}}}^{(a)} 〛 (ξ - ξ^{'}) \frac{\partial}{\partial v} {\bar{e}}_{z} (ξ^{'}, v) d ξ^{'} \\ = - i ω μ_{r} {\bar{h}}_{y} (ξ, v) \end{array}

where Ū(ξ), 〚ġ^(p)〛, and 〚ġ̄^(a)〛(ξ) are (2N + 1) × (2N + 1) square matrices whose (n,m)-entries are given by

{(\bar{U} (ξ))}_{n, m} = δ_{n, m} α_{n} (ξ)

{(〚 {\dot{g}}^{(p)} 〛)}_{n, m} = \frac{1}{d} \int_{- d / 2}^{d / 2} {\dot{g}}^{(p)} (u) e^{- i (n - m) k_{d} u} d u

{(〚 {\bar{\dot{g}}}^{(a)} 〛 (ξ))}_{n, m} = \frac{1}{d} \int_{- \infty}^{\infty} g^{(a)} (u) e^{- i α_{n - m} (ξ) u} d u .

The symbol δ_n,m stands for Kronecker’s delta. The integrands appeared in Eqs. (10) and (12) are products of pseudo-periodic functions with no concurrent jump discontinuities, and the Laurent rule [12] has been used to derive Eqs. (15) and (17).

Considering the periodicity in terms of ξ, Eqs. (15)–(17) have to be satisfied for arbitrary ξ ∈ (−k_d/2, k_d/2]. However, we take L sample points inside the Brillouin zone and assume that Eqs. (15)–(17) are satisfied at these sample points. Also the integrations in Eqs. (15) and (17) are approximated by an appropriate numerical integration scheme using the same sample points. Let ${ξ_{l}}_{l = 1}^{L}$ and ${w_{l}}_{l = 1}^{L}$ be the sample points and the weights chosen by a numerical integration scheme. Then we may obtain the following relations:

(i \tilde{U} - \tilde{G} \frac{d}{d v}) {\tilde{h}}_{y} (v) - \frac{d}{d v} {\tilde{h}}_{x} (v) = - i ω ε_{r} {\tilde{e}}_{z} (v)

\frac{d}{d v} {\tilde{e}}_{z} (v) = i ω μ_{r} {\tilde{h}}_{x} (v)

(i \tilde{U} - \tilde{G} \frac{d}{d v}) {\tilde{e}}_{z} (v) = - i ω μ_{r} {\tilde{h}}_{y} (v)

with

{\tilde{e}}_{z} (v) = (\begin{matrix} {\tilde{e}}_{z} (ξ_{1}, v) \\ ⋮ \\ {\tilde{e}}_{z} (ξ_{L}, v) \end{matrix})

\tilde{U} = (\begin{matrix} \bar{U} (ξ_{1}) & 0 \\ ⋱ \\ 0 & \bar{U} (ξ_{L}) \end{matrix})

\tilde{G} = (\begin{matrix} {\bar{G}}_{1, 1} & \dots & {\bar{G}}_{1, L} \\ ⋮ & ⋱ & ⋮ \\ {\bar{G}}_{L, 1} & \dots & {\bar{G}}_{L, L} \end{matrix})

{\bar{G}}_{l, l^{'}} = δ_{l, l^{'}} 〚 {\dot{g}}^{(p)} 〛 + \frac{w_{l^{'}}}{k_{d}} 〚 {\bar{\dot{g}}}^{(a)} 〛 (ξ_{l} - ξ_{l^{'}})

where the definition of the column matrices h̃_x(v) and h̃_y(v) are similar to ẽ_z(v). After a simple calculation, Eqs. (21)–(23) yield the following coupled differential-equation set:

(\begin{matrix} {\tilde{e}}_{z} (v) \\ - i \frac{d}{d v} {\tilde{e}}_{z} (v) \end{matrix}) = - i M_{r} \frac{d}{d v} (\begin{matrix} {\tilde{e}}_{z} (v) \\ - i \frac{d}{d v} {\tilde{e}}_{z} (v) \end{matrix})

with

M_{r} = (\begin{matrix} - {(k_{r}^{2} I - {\tilde{U}}^{2})}^{- 1} (\tilde{U} \tilde{G} + \tilde{G} \tilde{U}) & {(k_{r}^{2} I - {\tilde{U}}^{2})}^{- 1} ({\tilde{G}}^{2} + I) \\ I & 0 \end{matrix})

where I denotes the identity matrix and the superscript “−1” stands for the matrix inverse. The general solution to the coupled differential equation set (28) can be obtained by solving the eigenvalue-eigenvector problems because the matrix of coefficients M_r is constant. The 2L(2N + 1) eigenvalues can be divided into two sets, each containing L(2N + 1) eigenvalues. The first set contains the negative real eigenvalues and the complex eigenvalues with positive imaginary parts, and the second set contains those with the opposite signs. We denote the reciprocals of the eigenvalues of M_r by

{η_{r, n}}_{n = 1}^{2 L (2 N + 1)}

, in which

{η_{r, n}}_{n = 1}^{L (2 N + 1)}

correspond to the first set and

{η_{r, n}}_{n = L (2 N + 1) + 1}^{2 L (2 N + 1)}

correspond to the second set. Let p_r,n denote the eigenvector of M_r associating with the eigenvalue 1/η_r,n. Then the matrix for diagonalization is constructed as

(\begin{array}{l} {\tilde{P}}_{r, 11} & {\tilde{P}}_{r, 12} \\ {\tilde{P}}_{r, 21} & {\tilde{P}}_{r, 22} \end{array}) = (\begin{matrix} p_{r, 1} & \dots & p_{r, 2 L (2 N + 1)} \end{matrix})

and the general solution to the coupled differential-equation set (28) is written in the following form:

(\begin{matrix} {\tilde{e}}_{z} (v) \\ - i \frac{d}{d v} {\tilde{e}}_{z} (v) \end{matrix}) = (\begin{array}{l} {\tilde{P}}_{r, 11} & {\tilde{P}}_{r, 12} \\ {\tilde{P}}_{r, 21} & {\tilde{P}}_{r, 22} \end{array}) (\begin{matrix} {\tilde{a}}_{e, r}^{(-)} (v) \\ {\tilde{a}}_{e, r}^{(+)} (v) \end{matrix})

where the column matrices

{\tilde{a}}_{e, r}^{(-)} (v)

and

{\tilde{a}}_{e, r}^{(+)} (v)

are the amplitudes of the eigenmodes propagating in the negative and the positive v-directions, respectively. The relation between the modal amplitudes at v = v′ and v = v″ is given by

(\begin{matrix} {\tilde{a}}_{e, r}^{(-)} (v^{'}) \\ {\tilde{a}}_{e, r}^{(+)} (v^{'}) \end{matrix}) = V_{r} (v^{'} - v^{″}) (\begin{matrix} {\tilde{a}}_{e, r}^{(-)} (v^{″}) \\ {\tilde{a}}_{e, r}^{(+)} (v^{″}) \end{matrix})

with

{(V_{r} (v))}_{n, m} = δ_{n, m} e^{i η_{r, n} v} .

The covariant component of the magnetic field in terms of u is given by H_t (u, v) = H_x(u, v) + ġ(u)H_y(u, v), which gives the tangential component of the magnetic field on the grating surface v = 0. From Eqs. (22), (23), and (31), the generalized Fourier coefficients of Ē_z(u;ξ_l, v) and H̄_t(u;ξ_l, v) are expressed in the following form:

(\begin{matrix} {\tilde{e}}_{z} (v) \\ {\tilde{h}}_{t} (v) \end{matrix}) = (\begin{array}{l} {\tilde{P}}_{r, 11} & {\tilde{P}}_{r, 12} \\ {\tilde{Q}}_{e, r, 1} & {\tilde{Q}}_{e, r, 2} \end{array}) (\begin{matrix} {\tilde{a}}_{e, r}^{(-)} (v) \\ {\tilde{a}}_{e, r}^{(+)} (v) \end{matrix})

with

{\tilde{Q}}_{e, r, q} = - \frac{1}{ω μ_{r}} [\tilde{G} \tilde{U} {\tilde{P}}_{r, 1 q} - ({\tilde{G}}^{2} + I) {\tilde{P}}_{r, 2 q}]

for q = 1,2.

For the TM-polarization, following the same process, the generalized Fourier coefficients of H̄_z(u;ξ_l, v) and Ē_t (u;ξ_l, v) are expressed in the following form:

(\begin{matrix} {\tilde{h}}_{z} (v) \\ {\tilde{e}}_{t} (v) \end{matrix}) = (\begin{matrix} {\tilde{P}}_{r, 11} & {\tilde{P}}_{r, 12} \\ {\tilde{Q}}_{h, r, 1} & {\tilde{Q}}_{h, r, 2} \end{matrix}) (\begin{matrix} {\tilde{a}}_{h, r}^{(-)} (v) \\ {\tilde{a}}_{h, r}^{(+)} (v) \end{matrix})

with

{\tilde{Q}}_{h, r, q} = \frac{1}{ω ε_{r}} [\tilde{G} \tilde{U} {\tilde{P}}_{r, 1 q} - ({\tilde{G}}^{2} + I) {\tilde{P}}_{r, 2 q}]

for q = 1,2. It is worth noting that the coefficient matrix M_r is independent of the polarization and, therefore, P̃_r,pq for p, q = 1, 2 are also independent of the polarization.

The general solutions separately obtained in the regions s and c can be matched at the grating surface v = 0 by using the boundary conditions, which are given by the continuities of the tangential components of the fields. Then the scattering matrix that relates the amplitudes of the incoming and outgoing fields is derived as

(\begin{matrix} {\tilde{a}}_{f, s}^{(+)} (+ 0) \\ {\tilde{a}}_{f, c}^{(-)} (- 0) \end{matrix}) = (\begin{matrix} {\tilde{S}}_{f, 11} & {\tilde{S}}_{f, 12} \\ {\tilde{S}}_{f, 21} & {\tilde{S}}_{f, 22} \end{matrix}) (\begin{matrix} {\tilde{a}}_{f, s}^{(-)} (+ 0) \\ {\tilde{a}}_{f, c}^{(+)} (- 0) \end{matrix})

with

{\tilde{S}}_{f, 11} = - {({\tilde{Q}}_{f, c, 1} {\tilde{P}}_{c, 11}^{- 1} {\tilde{P}}_{s, 12} - {\tilde{Q}}_{f, s, 2})}^{- 1} ({\tilde{Q}}_{f, c, 1} {\tilde{P}}_{c, 11}^{- 1} {\tilde{P}}_{s, 11} - {\tilde{Q}}_{f, s, 1})

{\tilde{S}}_{f, 12} = {({\tilde{Q}}_{f, c, 1} {\tilde{P}}_{c, 11}^{- 1} {\tilde{P}}_{s, 12} - {\tilde{Q}}_{f, s, 2})}^{- 1} ({\tilde{Q}}_{f, c, 1} {\tilde{P}}_{c, 11}^{- 1} {\tilde{P}}_{c, 12} - {\tilde{Q}}_{f, c, 2})

{\tilde{S}}_{f, 21} = {\tilde{P}}_{c, 11}^{- 1} ({\tilde{P}}_{s, 12} {\tilde{S}}_{f, 11} + {\tilde{P}}_{s, 11})

{\tilde{S}}_{f, 22} = {\tilde{P}}_{c, 11}^{- 1} ({\tilde{P}}_{s, 12} {\tilde{S}}_{f, 12} - {\tilde{P}}_{c, 12})

for f = e,h. From Eq. (31), the amplitudes of the eigenmodes are written as

(\begin{matrix} {\tilde{a}}_{f, r}^{(-)} (v) \\ {\tilde{a}}_{f, r}^{(+)} (v) \end{matrix}) = {(\begin{array}{l} {\tilde{P}}_{r, 11} & {\tilde{P}}_{r, 12} \\ {\tilde{P}}_{r, 21} & {\tilde{P}}_{r, 22} \end{array})}^{- 1} (\begin{matrix} \tilde{ψ} (v) \\ - i \frac{d}{d v} \tilde{ψ} (v) \end{matrix})

where ψ̃(v) denotes ẽ_z(v) for the TE-polarization and h̃_z(v) for the TM-polarization. The column matrices ψ̃(v) and

- i \frac{d}{d v} \tilde{ψ} (v)

are generated the generalized Fourier coefficients of the transformed field, which can be written in

{\bar{ψ}}_{n} (ξ, v) = \frac{1}{d} \int_{- \infty}^{\infty} ψ (u, v + g (u)) e^{- i α_{n} (ξ) u} d u

- i \frac{\partial}{\partial v} {\bar{ψ}}_{n} (ξ, v) = - \frac{i}{d} \int_{- \infty}^{\infty} \frac{\partial}{\partial v} ψ (u, v + g (u)) e^{- i α_{n} (ξ) u} d u

for the original field ψ(x,y). The amplitudes of the incoming fields

{\tilde{a}}_{f, s}^{(-)} (+ 0)

and

{\tilde{a}}_{f, c}^{(+)} (- 0)

may be obtained by using these relations for known incident fields.

4. Numerical experiments

To validate the present formulation, we consider sinusoidal profile gratings that have local collapses of the structural perfection, and show the results of numerical experiments. The periodic part of the surface profile function g^(p)(x) is given by

g^{(p)} (x) = h_{1} + (h_{2} - h_{1}) {sin}^{2} (\frac{π}{d} x)

and the fields are excited by a line-source situated parallel to the z-axis at (x,y) = (x₀, y₀) for y₀ > h₂. The grating parameters are chosen as follows: d = 0.6λ₀, h₁ = −0.5d, h₂ = 0, ε_s = ε₀, ε_c = (1.3 + i7.6)² ε₀ (a typical value of Aluminum at λ₀ = 632.8nm), and μ_s = μ_c = μ₀.

4.1. Treatment of line-source excitation

Since the line-source is located in the region s, the incident field ψ⁽ⁱ⁾(x,y) is expressed as

ψ^{(i)} (x, y) = H_{0}^{(1)} (k_{s} \sqrt{{(x - x_{0})}^{2} + {(y - y_{0})}^{2}})

for y > g(x), where

H_{0}^{(1)}

denotes the zeroth-order Hankel function of the first kind. For the line-source, the direct integrations in Eqs. (44) and (45) are not easy, and we rewrite these equations as

\begin{array}{l} {\bar{ψ}}_{n}^{(i)} (ξ, v) & = & \frac{1}{d} \int_{- \infty}^{\infty} ψ^{(i)} (u, v + g^{(p)} (u)) e^{- i α_{n} (ξ) u} d u \\ + \frac{1}{d} \int_{a_{1}}^{a_{2}} (ψ^{(i)} (u, v + g (u)) - ψ^{(i)} (u, v + g^{(p)} (u))) e^{- i α_{n} (ξ) u} d u \end{array}

\begin{array}{l} - i \frac{\partial}{\partial v} {\bar{ψ}}_{n}^{(i)} (ξ, v) = - \frac{i}{d} \int_{- \infty}^{\infty} \frac{\partial}{\partial v} ψ^{(i)} (u, v + g^{(p)} (u)) e^{- i α_{n} (ξ) u} d u \\ - \frac{i}{d} \int_{a_{1}}^{a_{2}} \frac{\partial}{\partial v} (ψ^{(i)} (u, v + g (u)) - ψ^{(i)} (u, v + g^{(p)} (u))) e^{- i α_{n} (ξ) u} d u . \end{array}

Here, we use the Fourier integral representation of the Hankel function of the first kind of order zero [13]:

H_{0}^{(1)} (k_{s} ρ (x, y)) = \frac{1}{π} \int_{- \infty}^{\infty} \frac{1}{\sqrt{k_{s}^{2} - ζ^{2}}} e^{i (ζ x + \sqrt{k_{s}^{2} - ζ^{2}} | y |)} d ζ

and express the periodic function exp(iβ(h₂ − g^(p)(u))) in the Fourier series as

e^{i β (h_{2} - g^{(p)} (u))} = \sum_{m = - \infty}^{\infty} c_{m} (β) e^{i m k_{d} u} .

Then, the first integrals on the right-hand sides of Eqs. (48) and (49) can be written as

\begin{array}{l} \frac{1}{d} \int_{- \infty}^{\infty} ψ^{(i)} (u, v + g^{(p)} (u)) e^{- i α_{n} (ξ) u} d u \\ = \sum_{m = - \infty}^{\infty} c_{n - m} (β_{m} (ξ)) \frac{2}{d β_{m} (ξ)} e^{- i [α_{m} (ξ) x_{0} - β_{m} (ξ) (y_{0} - h_{2})]} e^{- i β_{m} (ξ) v} \end{array}

\begin{array}{l} - \frac{i}{d} \int_{- \infty}^{\infty} \frac{\partial}{\partial v} ψ^{(i)} (u, v + g^{(p)} (u)) e^{- i α_{n} (ξ) u} d u \\ = - \frac{2}{d} \sum_{m = - \infty}^{\infty} c_{n - m} (β_{m} (ξ)) e^{- i [α_{m} (ξ) x_{0} - β_{m} (ξ) (y_{0} - h_{2})]} e^{- i β_{m} (ξ) v} \end{array}

with

β_{m} (ξ) = \sqrt{k_{s}^{2} - α_{m} {(ξ)}^{2}}

near the corrugated surface. The periodic part of surface profile function is given by Eq. (46) and, then, the Fourier coefficients appeared in Eq. (51) are given by

c_{m} (β) = i^{m} J_{m} (β \frac{h_{1} - h_{1}}{2}) e^{i β \frac{h_{2} - h_{1}}{2}}

where J_m denotes the mth-order Bessel function. We truncate the infinite sums in Eqs. (52) and (53) by using the same truncation order N with Eq. (13) and obtain approximate values of the first integrals on the right-hand sides of Eqs. (48) and (49). On the other hand, the second integrals are with a finite interval, and can be evaluated by standard numerical integration.

4.2. Grating with a defect

The first example of the imperfectly periodic structure is the grating with a defect. Here, we define the aperiodic part of the surface profile function g^(a)(x) as

g^{(a)} (x) = {\begin{array}{l} (h_{2} - h_{1}) {cos}^{2} (\frac{π}{d} x) & for | x | \leq \frac{d}{2} \\ 0 & for | x | > \frac{d}{2} \end{array}

and, then, one groove of the original grating is filled by the substrate medium. We set the position of line-source at (x₀, y₀) = (d, 2d) and the observation point for convergence tests (x,y) = (0,d).

The obtained field intensities at the observation point are shown in Fig. 2. Figure 2(a) shows the convergence characteristics in terms of the number of sample points L. The truncation order for each generalized Fourier series expansion is set to N = 4 to calculate these results. The dotted curves are the results of the trapezoidal scheme that uses equidistant sample points and equal weights. Since the transformed field ψ̄(x;ξ,y) is not a smooth function of ξ, they converge very slowly. On the other hand, for the solid curves, we split the Brillouin zone at the Wood-Rayleigh anomalies, which is the non-smooth points, and the double exponential scheme [14, 15] is applied for the subintervals to decide the sample points and the weights. This scheme was used in Ref. [7], which presents RCWA with the help of PPFT for the lamellar gratings with defects, and showed significant improvement of convergence. The solid curves in Fig. 2(a) show that the same scheme is also effective to improve the convergence. The results of convergence test with respect to the truncation order N are shown in Fig. 2(b). The values are computed with L = 80, and the sample points and the weights are determined by applying the double exponential scheme for the subintervals. The convergence with respect to the truncation order N is very fast. Since the present formulation does not include the generalized Fourier series expansion of discontinuous function and the surface profile under consideration is smooth, the convergence speed is slightly faster than that shown in Fig. 2(b) of Ref. [7]. The distribution of the total field intensities near the defect are calculated with L = 80 and N = 4 by changing the observation point, and shown in Fig. 3. The positions of surfaces are indicated by the white dashed lines and the obtained results seem to be proper.

Fig. 2 Convergence of the field intensities at (x,y) = (0,d) for a sinusoidal profile grating with a defect.

Download Full Size | PDF

Fig. 3 Field intensities near a sinusoidal profile grating with a defect.

Download Full Size | PDF

We examine also the reciprocal property to validate the present formulation. We use the reciprocity error introduced in Ref. [7] here. Let ψ(x_p, y_p; x_q, y_q) denote the field observed at (x_p, y_p) for a line-source located at (x_q, y_q). Then, the reciprocity error for two points (x_p, y_p) and (x_q, y_q) is defined by

σ (x_{p}, y_{p}; x_{q}, y_{q}) = \frac{| ψ (x_{p}, y_{p}; x_{q}, y_{q}) - ψ (x_{q}, y_{q}; x_{p}, u_{p}) |}{| ψ (x_{p}, y_{p}; x_{q}, y_{q}) |} .

The reciprocity theorem requires that this function is zero when both (x_p, y_p) and (x_q, y_q) are located in the surrounding medium. We fix one point at (x_p, y_p) = (0, 2d) and the other point (x_q, y_q) is moved on the line y = d. The reciprocity errors at 101 equidistant points in −5d ≤ x ≤ 5d are calculated with L = 80 and N = 4 in the standard double-precision arithmetic, and the obtained results are shown in Fig. 4. The largest values are about 2.9 × 10⁻³ for the TE-polarization and 2.9 × 10⁻² for the TM-polarization. The largest value for the TM-polarization does not seem to be small enough but the related field intensity is very small (about 5.3×10⁻³) for this point. The second largest value for the TM-polarization is about 3.4 × 10⁻³, and it can be said that the reciprocity relation is well satisfied.

Fig. 4 Numerical results of the reciprocity test for a sinusoidal profile grating with a defect.

Download Full Size | PDF

4.3. Grating with period modulation

Next, we consider the grating in which the period is locally modulated. We define the aperiodic part of the surface profile function as

g^{(a)} (x) = {\begin{array}{l} - g^{(p)} (x) + g^{(p)} (\int_{0}^{x} \frac{d}{d + Δ (ζ)} d ζ) & for | x | \leq a \\ 0 & for | x | > a \end{array}

with

Δ (x) = \frac{d_{min} - d}{2} (1 + cos (\frac{π}{a} x)) .

The local period of the grating is modulated to be d + Δ(x) in the region −a ≤ x ≤ a. The local period deviation Δ(x) is a raised-cosine function and the minimum value of the local period is denoted by d_min. The parameter a is chosen to be a = d/[(d/d_min)^1/2 − 1], and the surface profile function g(x) is smooth even at x = ±a. Here, the positions of the line-source and the observation point are respectively chosen at (x₀, y₀) = (d, 2d) and (x, y) = (0, d), and the minimum local period is d_min = 0.5d.

Figure 5 shows the results of the same convergence tests as in Fig. 2 though the values in Fig. 5(a) are computed with N = 8. The convergence characteristics in terms of the sample point number L are similar to those for the grating with a defect. In contrast, the convergence characteristics in terms of the truncation order N are slower than that for the grating with a defect. Since the surface profile function g(x) under consideration has larger components in high spatial frequency range than that for the grating with a defect, the generalized Fourier series expression requires more terms to construct the original function accurately. The field intensity distributions computed with L = 80 and N = 8 and shown in Fig. 6 and the obtained results seem to be proper for both of the TE- and the TM-polarized cases. Figure 7 shows the results of the same reciprocity test with Fig. 2 though the values are computed with L = 80 and N = 8. The largest values are about 1.5 × 10⁻³ for the TE-polarization and 1.2 × 10⁻³ for the TM-polarization, and the reciprocity relation is well satisfied.

Fig. 5 Convergence of the field intensities at (x,y) = (0,d) for a sinusoidal profile grating with period modulation.

Download Full Size | PDF

Fig. 6 Field intensities near a sinusoidal profile grating with period modulation for the line-source excitation.

Download Full Size | PDF

Fig. 7 Numerical results of the reciprocity test for a sinusoidal profile grating with period modulation.

Download Full Size | PDF

5. Concluding remarks

This paper has dealt with the two-dimensional electromagnetic scattering from periodically corrugated surface in which structural periodicity is locally collapsed, and presented a spectral-domain formulation based on the C-method with the help of PPFT. The main problem on the spectral-domain analysis is generally summarized to the discretization scheme on the wavenumber space. In the present formulation, the discretization scheme can be considered only inside the Brillouin zone owing to the PPFT. To validate the present formulation, we considered two specific examples of imperfectly periodic structures: a sinusoidal profile grating with a defect and the same grating with period modulation, and showed some numerical results for a line-source excitation. A simple discretization scheme presented in Ref. [7] showed good convergences, and the convergence speeds in terms of the sample point numbers were similar for both examples. By contrast, the convergence characteristics in terms of the truncation order depend on the spatial frequency of the structure. The present formulation was also validated by the reciprocity tests.

This work has been partially supported by the Grant Agency of the Czech Republic (# P205/11/2137), by Structural Funds of the European Union and state budget of the Czech Republic ( NANOBASE project # CZ.1.07/2.3.00/20.0074), and by the IT4Innovations Centre of Excellence project, reg. no. CZ.1.05/1.1.00/02.0070.

References and links

1. T. Oonishi, T. Konishi, and K. Itoh, “Fabrication of phase only binary blazed grating with subwavelength structures designed by deterministic method based on electromagnetic analysis,” Jpn. J. Appl. Phys. 46, 5435–5440 (2007). [CrossRef]

2. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Modeling the Flow of Light (Princeton Univ. Press, Princeton, 1995).

3. C. Yang, K. Shi, P. Edwards, and Z. Liu, “Demonstration of a PDMS based hybrid grating and Fresnel lens (G-Fresnel) device,” Opt. Express 18, 23529–23534 (2010). [CrossRef] [PubMed]

4. J. Chandezon, M. T. Dupuis, G. Cornet, and D. Maystre, “Multicoated gratings: a differential formalism applicable in the entire optical region,” J. Opt. Soc. Am. 72, 839–846 (1982). [CrossRef]

5. R. Petit, ed., Electromagnetic Theory of Gratings (Springer-Verlag, Berlin, 1980). [CrossRef]

6. K. Watanabe and K. Yasumoto, “Two-dimensional electromagnetic scattering of non-plane incident waves by periodic structures,” Prog. Electromagnetic Res. PIER 74, 241–271 (2007). [CrossRef]

7. K. Watanabe, J. Pištora, and Y. Nakatake, “Rigorous coupled-wave analysis of electromagnetic scattering from lamellar grating with defects,” Opt. Express 19, 25799–25811 (2011). [CrossRef]

8. K. Knop, “Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves,” J. Opt. Soc. Am. 68, 1206–1210 (1978). [CrossRef]

9. M. G. Moharam and T. K. Gaylord, “Diffraction analysis of dielectric surface-relief gratings,” J. Opt. Soc. Am. 72, 1385–1392 (1982). [CrossRef]

10. E. Popov, M. Nevière, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A 19, 33–42 (2002). [CrossRef]

11. K. Watanabe, “Numerical integration schemes used on the differential theory for anisotropic gratings,” J. Opt. Soc. Am. A 19, 2245–2252 (2002). [CrossRef]

12. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]

13. W. C. Chew, Waves and Fields in Inhomogeneous Media (Van Nostrand Reinhold, New York, 1990).

14. H. Takahasi and M. Mori, “Double exponential formulas for numerical integration,” Publ. RIMS, Kyoto Univ. 9, 721–741 (1974). [CrossRef]

15. P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, 2nd ed. (Academic Press, New York, 1984).

Coordinate transformation formulation of electromagnetic scattering from imperfectly periodic surfaces

Abstract

1. Introduction

2. Settings of the problem

3. Formulation

4. Numerical experiments

4.1. Treatment of line-source excitation

4.2. Grating with a defect

4.3. Grating with period modulation

5. Concluding remarks

References and links

Cited By

Figures (7)

Equations (59)

Optics Express