1. Introduction

Fourier-expansion-based electromagnetic theories treating the optical response of periodic structures such as diffraction gratings (hereafter referred to as the coupled wave theory [1, 2, 3]) or photonic crystals [4, 5] have been the subject of considerable development for past several decades. In the case of discontinuous structures the method suffered from poor convergence until the correct Fourier factorization (FF) rules were introduced by Li [6]. Those rules have then been successfully applied to anisotropic [7], slanted [8], arbitrary-relief [9], two-dimensional (2D) gratings [10], and their various combinations [11, 12, 13], and are also used for other applications such as nonperiodic cylindrical devices [14, 15] or nonlinear optical systems [16].

As regards 2D structures [17, 18, 19], considerable troubles appeared in the case of elements with circular (or elliptical) cross-sections for which Li’s method of successive one-dimensional (1D) expansions [10] was not convenient. For this purpose David et al. [20] applied a method treating independently the tangential and normal components of fields on the edges of circular (or elliptical) holes, forming 2D photonic crystals. Similarly, Schuster and colleagues [21, 22] employed the normal vector method to 2D gratings with various distributions of polarization bases. However, these approaches, always dealing with linear polarizations, ignored the fact that the transformation matrix between the Cartesian and the normal/tangential component bases of polarization became discontinuous at the center and on the boundaries of the periodic cell, which slowed down the convergence of the numerical implementation.

In this article the theory is reformulated via FF with complex polarization bases, which avoids the points of discontinuity of the polarization transformation by employing generally elliptical polarization bases, simply by using complex-valued Jones matrices [23]. The FF approach is incorporated into the Airy-like series propagation algorithm (which can be regarded as a physically transparent variant of the scattering-matrix algorithm [24, 25, 26]), following a 1D analogy described in [27]; the methodology for obtaining matrix elements is in [18].

 

Fig. 1. Configuration of the optical problem (a) and propagation algorithm (b).

Download Full Size | PDF

2. Theoretical models

The coupled wave theory of 2D gratings is based on treating the equations (note that the space coordinates and the magnetic field [H_x; H_y] are scaled)

in the Fourier space by expanding the bi-periodic relative permittivity function ε(x,y) and the transverse (to the z axis) components of quasi-periodic electric field E _x(y)

(and analogously H_x, H_y), where ∂_x = ∂/∂_x etc., ∇_t = [∂_x,∂_y], p_m = sinϑⁱcosφⁱ + mλ/Λ, q_n = sinϑⁱsinφⁱ + nλ/Λ, with λ denoting the wavelength, Λ the grating period, ϑⁱ and φⁱ the spherical angles of the incident wave, and s the eigennumber of the wave equation (1). The configuration of the optical problem is displayed in Fig. 1(a), and the principle of the Airy-like series propagation algorithm is depicted in Fig. 1(b), where the propagation matrix P is determined from the eigensolutions of Eq. (1) and where the reflection and transmission matrices on interfaces R _jk, T _jk follow from boundary-matching conditions using Eq. (2).

In the Fourier space the operators ∂_x, ∂_y become diagonal matrices -i p, -i q with particular arrangements of the p_m, q_n values on their diagonals [18], whereas the multiplication εE_x (or εE_y) can be treated by either the Laurent factorization rule [εE_x] = [[ε]][E_x] (matrix version of the convolution rule between the Fourier coefficients of the multiplied functions) or by the inverse rule [εE_x ] = [[1/ε]]^-1 [E_x], where [f] denotes a column vector of the Fourier coefficients of a function f, while [[g]] denotes a matrix composed of the Fourier coefficients of a function g [18]. In this article we will define three models according to applying different factorization rules; first (Model A) assumes the Laurent rules [εE_x ] = [[ε]] [E_x] and [εE_y] = [[ε]] [E_y].

According to [6] the Laurent rule is valid (for the reason of uniform convergence) for multiplied functions possessing no concurrent discontinuities, whereas the inverse rule is used for functions whose product is continuous. Obviously, neither the Laurent rule nor the inverse rule is correct for both products D_x = εE_x and D_y = εE_y, because both pairs of functions have concurrent discontinuities and both products D32 and D_y (the electrical displacement) are discontinuous as well. On the other hand, by a clever change of the polarization bases at all points (using a space-dependent Jones matrix transform F), we can treat independently the normal and tangential components of the fields by the correct rules, i.e., [D_u] = [[1/ε]]^-1[E_u] and [D_v] = [[ε]][E_v], where E_u, D_u are field components normal to the discontinuities of the permittivity function, while E_v, D_v are tangential; the reason for this is because E_v and D_u are continuous. This idea was applied in [9] for 1D gratings, and proposed in [2] and applied in [21] for 2D periodicity.

 

Fig. 2. Distribution of u in the periodic cell, corresponding to Models B (a) and C (b), and the radial dependence of the ellipticity along two lines ϕ = 0 and ϕ = 30° (c).

Download Full Size | PDF

According to [21], the simplest distribution of the matrix F, whose columns are the basis polarization vectors u = [ξ,; ζ], v = [-ζ ^*; ξ ^*] (mutually orthogonal [23]), can be defined as the rotation by a polar angle ϕ, i.e., ξ = cos ϕ, ζ = sinϕ, within a single cell where we define polar coordinates re^iϕ = x + i_y, and then by periodic repeating over the entire 2D space. The distribution of the basis vector u in the periodic cell is depicted in Fig. 2(a), from where it is obvious that the matrix function F(x,y) has no concurrent discontinuities with the electric field, so that we can use the Laurent rule for the transformation of polarization

which yields the 2D factorization rule directly applicable to Eqs. (1) and (2),

Unfortunately, this approach (here referred to as Model B) only deals with linear polarizations and thus suffers from the fact that the matrix function F(x,y) has a singularity at the center of the periodic cell and other discontinuities on the cell boundaries, which slows down the convergence of the numerical implementation, as will be evidenced below.

To avoid the discontinuities, we develop the method of Fourier factorization with complex polarization bases (referred to as Model C) by using generally elliptic polarizations, so that the matrix function F(x,y) becomes completely continuous. We still use Eqs. (5) and (6) but with

where

is the space distribution (in the polar coordinates defined in one periodic cell) of the azimuth and the ellipticity of the polarization ellipse of the basis vector u, as depicted in Fig. 2(b), with R denoting the radius of the cylindrical element and D(ϕ)= Λ/2max(∣cos ϕ∣, ∣sinϕ∣) being the distance from the cell’s center to its edge. In Eq. (7) the Jones vector on the right represents a polarization ellipse (with ellipticity η) oriented along the x coordinate, the matrix in the middle rotates this polarization by the azimuth θ, and the factor e^iθ preserves the continuity of the phase at the center and on the boundaries of the cell. Similarly to Model B, the azimuth of the polarization ellipse is constant on the lines coming from the cell’s center, but now the ellipticity is zero (corresponding to a linear polarization) only on the dot’s edges, has the maximum value (π/4 for the circular polarization) at the cell’s center and on its boundaries, and is continuously varying (with a smooth sine dependence) in the intermediate ranges [Fig. 2(c)]. The continuity at the center and boundaries of the periodic cell can be easily checked by evaluating lim_r→0 u = lim_r→D(ϕ) u = (1/2)^1/2[1; i], which is independent of ϕ. Hence we obtain a smooth and completely continuous matrix function F(x,y) corresponding to Model C.

3. Numerical example and discussion

Let us examine the numerical performance of this approach for a 2D grating made as bi-periodically arranged holes patterned on the top of a quartz substrate. Suppose cylindrical holes of the diameter 2R = 200 nm, depth d = 100 nm, square periodicity Λ = 300 nm, an incident plane wave with the wavelength λ = 500 nm, the angle of incidence ϑⁱ = 60°, and the angle of the plane of incidence φⁱ = 0 [the configuration of Fig. 1(a)], with the values of permittivity ε _a = 1 for vacuum inside holes and ε _b = 2.138 for quartz [28].

In Fig. 3 we compare calculations using Models A (dashed curves with crosses), B (dotted curves with empty squares), and C (solid curves with filled circles). We plot the energy reflectance in the specular reflection [Figs. 3(a) and 3(b)] and in a first-diffracted order [Figs. 3(c) and 3(d)], more precisely the diffraction efficiencies in the [0,0] and [-1,0] orders of reflection. The efficiencies are determined for the incident s and p polarizations, denoted R_ss ^(m,n) and R_pp ^(m,n) for the [m,n] order. Note that in the configuration of Fig. 1(a) the mixed-polarization (sp and ps) efficiencies of the two diffracted orders of interest are zero. We display the dependences of the quantities according to the maximum indices M = max ∣m∣ and N = max ∣n∣ of the Fourier harmonics retained inside the periodic medium. For simplicity we keep M = N, so that the order of the [[ε]]-like matrices is (2N + 1)². Due to the limited memory of the computer used, calculations were possible to perform for N ≤ 16.

As obvious from the plotted curves, all the three models converge to the same limit, but with considerably different convergence performances. While the values yielded by Model B possess precision of about 10^-4 for N ≥ 7, Model C yields precision of about 10^-5 for the same N (except the sensitive case of R_pp ^(0,0) with particularly low reflection). The values yielded by Model A are even one order worse than those of Model B, so that they are not visible in Figs. 3(a)–3(c); therefore, insets are included with longer ranges for comparison.

4. Conclusion

The coupled wave theory of bi-periodic diffraction gratings was reformulated by employing the FF method with complex polarization bases, which enabled considerable enhancement of the numerical capabilities, as evidenced from the comparison with previous approaches. For most practical purposes (e.g., for studying the tendency of the optical quantities while some parameter is varied, or for comparison with an experiment affected by measurement errors) the new method provides sufficient precision for N = 4 (in the presented configuration), which considerably saves the computer memory and the time of calculations. It is worth pointing out that the essential difference between the presented Model C and the previous Model B is (from the mathematical viewpoint) the fact that the matrix transformation of polarization in the former method contains complex-valued elements, so that this generalization is surprisingly simple.

 

Fig. 3. Convergence properties of diffraction efficiencies of periodically arranged cylindrical holes in the [0,0] (a, b) and [-1,0] (c, d) orders of diffraction for the s (a, c) and p (b, d) polarizations.

Download Full Size | PDF

The method can also be advantageously applied to photonic band calculations [by reformulating the permittivity matrix in the formula ε ^-1∇ × (∇ × E) = (ω/c)² E], nonperiodic cylindrical devices (by using a large periodic cell), etc., and can be straightforwardly generalized to elements with arbitrary cross-sections [by replacing the constant value of the radius R in Eq. (9) with a function R(ε) and by generalizing the azimuth θ(r,ε) so that the linear polarization on the edge of the element becomes again normal to it].