## Abstract

We describe an efficient implementation of the finite-difference time-domain (FDTD) method as applied to lightwave propagation through periodic media with arbitrary anisotropy (birefringence). A permittivity tensor with non-diagonal elements is successfully integrated here with periodic boundary conditions, bounded computation space, and the split-field update technique. This enables modeling of periodic structures using only one period even with obliquely incident light in combination with both monochromatic (sinusoidal) and wideband (time-domain pulse) sources. Comparisons with results from other techniques in four validation cases are presented and excellent agreement is obtained. Our implementation is freely available on the Web.

©2006 Optical Society of America

## 1. Introduction

Optical elements with anisotropic dielectric properties, such as waveplates and liquid crystal light valves, have achieved substantial practical use in recent years because of their ability to control the polarization of lightwaves. One of the most powerful approaches for the analysis of planar anisotropic structures is the 2×2 Jones method [1], which is widely used for most simple cases at normal incidence due to its concise expressions. The extended Jones method incorporates oblique incidence [2], and a more widely applicable method was suggested by Berreman [3]. These matrix-type solvers, however, are limited to one-dimensional structures since all are based on a stratified medium approach. In this paper, we describe a more flexible and general approach based on the finite-difference time-domain (FDTD) algorithm for the analysis of arbitrary, anisotropic periodic structures at oblique incidence.

Since the early work by Yee [4], FDTD solvers have been developed as one of the most effective tools for computational electromagnetics problems. FDTD methods present a robust and powerful approach to directly solve Maxwell’s curl equations both in time and space. In particular, FDTD techniques have an advantage for geometrically complex elements over frequency domain methods, such as the method of moments (MoM) [5], because of its ability to define arbitrary shapes and inhomogeneous properties on the grid-space. Another advantage of the FDTD techniques is its capability to visualize the real-time pictures of the electromagnetic wave.

Consider a two-dimensional, periodic structure in anisotropic media as shown in Fig. 1(a). When the structure is periodic, the problem dimension can be reduced to one period by applying periodic boundary conditions. Even though the structure is geometrically two dimensional, all *x*, *y*, and *z* components of each field variable should be considered because of the anisotropy of the media, which can produce coupling between all three field components. We employ a bi-dimensional Yee grid [6], which is a projection of the original Yee grid in 3D, on the *x*-*z* plane as illustrated in Fig. 1(b).

While the FDTD solution is straightforward for the normal incidence case, a difficulty arises for the obliquely incident case; phase variation across a unit period of the structure. This phase difference leads in time-advance or delay, which cannot be solved directly using time-domain techniques. There are a number of ways to solve for periodic boundary conditions at oblique incidence [7, 8]. The split-field update method [9] is one of the most successful methods to circumvent these phase variations using field transformation techniques. However, the previous work was limited to material properties with diagonal tensors. We have extended the split-field FDTD method for anisotropic materials.

Another important feature in FDTD methods is a proper absorbing boundary condition to truncate the simulation space without artificial reflections. The ‘uniaxial perfectly matched layer’ (UPML) can be implemented in the split-field technique as absorbing boundaries [10].

In this paper, we apply the extended split-field FDTD method to analyze periodic, anisotropic optical elements at a general angle of incidence. We limit ourselves to a 2D structure, but we note that this method can be extended to 3D in a straightforward manner. Section 2 describes the details of FDTD implementation for periodic anisotropic structures. In addition, we present a fully-vectorial near-field to far-field transformation to improve the accuracy for the highly off-axis propagation modes. Section 3 demonstrates the developed FDTD modeling of four different optical elements of anisotropic materials. Comparisons of the FDTD results with other methods, if applicable, are presented. Conclusions are drawn in Section 4.

## 2. FDTD implementation for periodic anisotropic structures

Maxwell’s equations for a nonconducting anisotropic medium can be written in the phasor form as follows:

where *ω* is the angular frequency, *ε*
_{0} and *µ*
_{0} are the permittivity and the permeability of free space, respectively, and the tilde denotes a tensor. Materials in this paper are assumed nonmagnetic so *$\tilde{\mu}$*=**I** (**I** is a 3×3 unity matrix).

If the material has linear dielectric properties, only three dielectric constants (*ε*
_{1,2,3}) and three
geometric angles (*α*,*β*, *γ*) are necessary to specify the full tensor description of *$\tilde{\epsilon}$*:

where *ε*
_{1,2,3} are the relative dielectric constants corresponding to the *principal axes* (**e**
_{1},**e**
_{2},**e**
_{3}) of the coordinate system (called the *principal system*, where the permittivity tensor is diagonal) [11]. **T** is the transformation matrix given by

where *α*, *β*, and *γ* are Euler anlges [12]. It is often convenient to introduce the impermittivity tensor *$\tilde{\kappa}$*=*$\tilde{\epsilon}$*
^{-1} and Eq. (1) can be rewritten as

#### 2.1. Split-field update method for periodic anisotropic media

Consider an incident planewave on a periodic structure at an angle *θ*_{inc}
from the *z*-axis as shown in Fig. 1(a). The phase difference between two boundaries at *x*=0 and *x*=Λ, where Λ is the period of the structure, can be accounted by introducing a new set of field variables:

where *c* is the speed of light for free space and *k*_{x}
=(*ω*/*c*) sin *θ*_{inc}
. Substituting *P*_{x,y,z}
and *Q*_{x,y,z}
into Eq. (4) and separating them into each field component yields:

where *κ*_{ab}
is the (*a,b)* element of *$\tilde{\kappa}$*, and

Once mapping the **P** and **Q** field components on the FDTD grid-space, *P*_{x,y,z}
and *Q*_{x,y,z}
have the same cell-to-cell field relationships as *E*_{x,y,z}
and *H*_{x,y,z}
at normal incidence. However, the off-diagonal elements of *$\tilde{\kappa}$* lead to additional time-derivative terms on the right-hand side of Eq. (6), which need a special treatment.

Similar to Ref. [9], the additional time-derivative terms can be eliminated by splitting field variables:

Substituting Eq. (8) into the left-hand side of Eq. (6) and (7) results in equations for “*a*” fields:

and

Unlike the normal FDTD update, the **P** and **Q** field variables must be updated simultaneously. As a result, storing all field values at every integer (*n*) and half-integer (
$n+\frac{1}{2}$
) time-step using separate variables is required. As an indicative example, the discrete expression for Eq. (9a) can be written as

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}+S{\kappa}_{\mathit{xy}}\left[({Q}_{x}{|}_{i+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}-{Q}_{x}{|}_{i+\frac{1}{2},k-\frac{1}{2}}^{n+\frac{1}{2}})-\left({Q}_{z}{|}_{i+1,k}^{n+\frac{1}{2}}{-Q}_{z}{|}_{i,k}^{n+\frac{1}{2}}\right)\right]$$

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}+S{\kappa}_{\mathit{xz}}({Q}_{y}{|}_{i+1,k}^{n+\frac{1}{2}}-{Q}_{y}{|}_{i,k}^{n+\frac{1}{2}})$$

where *S*=*c*Δ*t*/Δ*u* (Δ*u* and Δ*t* are the grid-spacing and the time-step, respectively). Note that we employ a square grid structure where Δ*x*=Δ*z*=Δ*u* for simplicity. The values of
${Q}_{x}{|}_{i+\frac{1}{2},k\pm \frac{1}{2}}^{n+\frac{1}{2}}$
, which are not available directly from the FDTD grid, must be interpolated [13] from neighboring grids as
$\frac{1}{2}({Q}_{x}{|}_{i+\frac{1}{2},k}^{n+\frac{1}{2}}+{Q}_{x}{|}_{i+\frac{1}{2},k\pm 1}^{n+\frac{1}{2}})$
. Similarly,
${Q}_{y,z}{|}_{i,k}^{n+\frac{1}{2}}$
and
${Q}_{y,z}{|}_{i+1,k}^{n+\frac{1}{2}}$
can be estimated;
${Q}_{y}{|}_{i,k}^{n+\frac{1}{2}}$
and
${Q}_{y}{|}_{i+1,k}^{n+\frac{1}{2}}$
require four neighboring quantities. These additional approximations still hold the second order of accuracy.

Once the “*a*” portions of the field variables are known, we can calculate the total fields. However, the remaining parts on the right-hand side of Eq. (8) still have temporal dependencies on the total fields, which must be removed properly. After some mathematical efforts, the total fields are found to be expressed using only the “*a*” fields as follows:

where *A*=1 - *κ*_{zz}
sin^{2}
*θ*_{inc}
, *B*=*κ*_{yz}
sin^{3} θ
_{inc}
, *C*=(1/*A*)*κ*_{yz}
sin^{2}
*θ*_{inc}
, and *D*=1 - (*B*/*A*)*κ*_{zy}
sin *θ*_{inc}
-*κ*_{yy}
sin^{2}
*θ*_{inc}
. The discrete version of Eq. (12a) can be found:

$$+\frac{\mathrm{sin}{\theta}_{\mathit{inc}}}{2D}({P}_{\mathit{ya}}{|}_{i,k}^{n+1}+{P}_{\mathit{za}}{|}_{i+1,k}^{n+1})$$

$$+\frac{B}{4\mathit{AD}}({P}_{\mathit{za}}{|}_{i,k-\frac{1}{2}}^{n+1}+{P}_{\mathit{za}}{|}_{i+1,k-\frac{1}{2}}^{n+1}+{P}_{\mathit{za}}{|}_{i,k+\frac{1}{2}}^{n+1}+{P}_{\mathit{za}}{|}_{i+1,k+\frac{1}{2}}^{n+1})$$

$$+\frac{C}{2D}({Q}_{\mathit{ya}}{|}_{i+\frac{1}{2},k-\frac{1}{2}}^{n+1}+{Q}_{\mathit{ya}}{|}_{i+\frac{1}{2},k+\frac{1}{2}}^{n+1}).$$

Applying the periodic boundary condition is now simply done by forcing field values at grid locations
$\left(i=\frac{\Lambda}{\Delta u}+\frac{1}{2}\right)$
and
$\left(i=\frac{1}{2}\right)$
to be the same as at
$\left(i=\frac{3}{2}\right)$
and
$\left(i=\frac{\Lambda}{\Delta u}-\frac{1}{2}\right)$
, respectively. We assume
$\frac{\Lambda}{\Delta u}$
is an integer for our convenience. The other two boundaries along the *x*-direction are modeled using UPML, which can be applied directly to the split-field method, as clearly implemented in Ref. [10].

#### 2.2. Numerical stability and dispersion

The numerical stability and dispersion are two main considerations of FDTD modeling with limited computation power. Since the FDTD method is based on sampling and updating the electric and magnetic field in time and space, wave propagation in the discrete grid-space may differ from it in the continuous space.

For a numerically stable FDTD method, the time step Δ*t* must be bounded to a special limit, often called the Courant factor (*S*=*c*Δ*t*/Δ*u*). A simple expression for the stability limit of the 2D split-field update method in free space was derived as a function of the incident angle *θ*_{inc}
[9]:

The exact stability limit for general anisotropic media can be found by solving an eleventh order polynomial [8], which is intractable since it also involves ten tightly-related model-dependent variables (including a permittivity tensor and incident angle). Instead of solving these complex equations, the same upper-bound of stability limit of Eq. (14) is still applicable to problems with anisotropic dielectric materials because of the following two reasons. First, the phase velocity in a non-dispersive dielectric medium is always slower than that in free space. Therefore, a system of dielectric media has a more relaxed stability condition than free space. Second, because we are only considering a 2D spatial grid and the additional spatial derivatives resulting from the non-diagonal permittivity tensor are implemented by the usual leap-frog scheme (with centered interpolations), the stability limit for the original 2D split-field method can be considered the upper bound for our case. All of our simulations, including those in Section 3, support these assertions.

The discrete grid-space and central difference approximations in the FDTD algorithm cause non-physical dispersion of the simulated waves, called numerical dispersion. For example, a spectral shift to longer wavelengths may appear as consequence of the numerical dispersion. Discussions on the numerical dispersion of the split-field techniques for homogeneous, isotropic materials can be found in Ref. [8]. However, the numerical dispersion relation for non-homogeneous media is not trivial, and depends on problem structures, grid-spacing, and time-step [14]. Obviously, one simple solution to reduce the impact of the numerical dispersion is smaller grid-spacing.

#### 2.3. Wideband source and vectorial far-field transformation

Two methods which create an input excitation are a continuous wave (CW) and a time-limited pulsed wave. The CW FDTD solution can be found in the steady state for sinusoidal illumination. A general form of the input sine wave is given by:

where **P**_{0}=**x̂**
*P*_{x}
exp(*jφ*_{x}
)+ **ŷ**
*P*_{y}
exp(*jφ*_{y}
). The complex amplitudes (*P*_{x}
and *P*_{y}
) and the phase difference (*φ*_{y}
-*φ*_{x}
) determine the polarization of the input wave. Note that the angle of incidence is already included in **P**_{0}, and therefore its propagation vector is always parallel to the *z*-axis. A pulsed planewave excitation can be applied to obtain a spectral response from a single FDTD simulation while the CW method requires an individual FDTD simulation for every frequency of interest. A gaussian pulse is widely used for pulsed FDTD [15, 16]:

where *ω*
_{0} is the peak-frequency, *t*
_{0} is a time-delay required to generate a smooth two-sided pulse, and *T*_{pulse}
determines the width of the pulse in the time-domain. The frequency information can be extracted by applying the discrete Fourier transformation during time-marching in the FDTD simulation.

The far-field information is often of interest in problems that involve periodicity. An infinitely periodic structure yields discrete diffracted orders (Floquet modes) in the region far from the structure. The far-field modes are governed by the so-called grating equation:

where *λ* is the wavelength of interest, Λ is the period, and *m* is an integer (often called Floquet modes). The far-field components can be calculated from the near-fields, *P*
_{(x,y,z),near}, by
applying the vectorial near-field to far-field transformation:

where ${E}_{\mathit{\text{TE}}\mathit{,}\mathit{\text{far}}}^{m}$
and ${E}_{\mathit{\text{TM}}\mathit{,}\mathit{\text{far}}}^{m}$
are the field components for the *TE* and *TM* modes of the *m*^{th}
order and *θ*′_{m}=sin^{-1}(*mλ*/Λ). We note that Eq. (18) is a more complete form of the field transformation over the conventional expression defined as

which is only valid for small angles and assumes a scalar description of the electric field.

## 3. FDTD validation

We apply the extended split-field FDTD method to analyze the optical properties of stratified and periodic anisotropic media at a general angle of incidence. Our discussion includes four different structures: a twisted-nematic liquid crystal cell, a planar slab of optically active media, a thin phase grating with rectangular grooves, and a special anisotropic grating known as a polarization grating (PG). To extract the far-field information, we apply our vectorial near-field to far-field transformation given in Eq. (18). FDTD results are compared with well-known analytical solutions or other numerical results such as the rigorous coupled-wave (RCW) theory and the Berreman method. Numerical parameters in Table 1 apply to all simulations in this paper. Note that *n*_{max}
is the largest index of refraction in the media and *λ*
_{0} is either the wavelength for monochromatic light or the peak-wavelength of a Gaussian pulse corresponding to *ω*
_{0}.

#### 3.1. Twisted-nematic liquid crystal cells

Consider a twisted-nematic (TN) liquid crystal (LC) layer between ideal cross polarizers as shown in Fig. 2. Adiabatic following occurs as light propagates through the LC domain when the nematic director manifests a slow twist (when *ϕ*_{twist}
≪πΔ*n*_{l}*d*/*λ*, where Δ*n*_{l}
is the LC birefringence, *d* is the LC thickness, and *ϕ*_{twist}
is the twist angle), referred to as the *Mauguin limit* [17]. The transmittance of the TN LC layer when placed between crossed polarizers was derived by Gooch and Tarry [18]:

where *ξ*_{l}
=Δ*n*_{l}*d*/*λ*. When the twist angle is 90°, the maximum transmittance occurs when:

The Berreman 4×4 matrix method [3] can be applied to analyze a TN LC layer, and will be used to verify the FDTD simulations. When a medium has one dimensional dependency and it can be easily stratified into a number of uniform layers, the Berreman method can yield accurate results.

A 90°-TN LC structure can be easily added in the FDTD grid-space by defining the permittivity tensor of Eq. (4) with *ε*
_{1,3}=${n}_{\perp}^{2}$, *ε*
_{2}=${n}_{\parallel}^{2}$, *α*(*z*)=(*z*/*d*)*ϕ*_{twist}
, and *β*=*γ*=0 for in-plane orientation of the LC director as

The modeled LC material parameters are *n*
_{⊥}=1.5,*n*
_{∥}=1.7, Δ*n*_{l}
=*n*
_{∥}-*n*
_{⊥}=0.2, and *ϕ*_{twist}
=90°. We assumed ideal/symbolic polarizers instead of implementing them in the FDTD grid; we excited a linearly polarized source pulse with half the incident intensity at the position of the first polarizer and we collected the near-field values of a single field component after the second (analyzing) polarizer. In addition, gradient-index anti-reflection (AR) coatings [19] are applied to minimize the effect of air-LC interfaces. These AR layers are implemented at both planar air-LC boundaries with a variable index distribution as

where *n*
_{1} is the refractive index of the incoming media, *n*
_{2} is the index for the outgoing media, and *t* is the AR thickness; we set *t*=*λ*
_{0} for the FDTD simulations. When the medium is anisotropic, we used the average index of refraction.

To analyze the optical properties of the TN LC cell, we calculated the Stokes parameters of light immediately after the LC layer and before the second polarizer. Fig. 3(a) shows the calculated polarization state of outgoing light from the TN LC cell in terms of the normalized Stokes parameters, *S*′_{1,2,3}=*S*
_{1,2,3}/*S*
_{0} where *S*
_{0} is the light intensity. Only the field component parallel to the analyzing polarizer can be transmitted through the analyzer. For our case, a high transmission can occur when *S*′_{1}≃-1. The transmittance through the TN LC cell is presented in Fig. 3(b). FDTD results are normalized by an AR-coated dielectric slab with the average index of refraction of the TN LC cell. An excellent agreement between FDTD and Berreman methods is found. To quantify numerical error, we calculated the mean-square error. For the transmittance of the TN LC cell, the mean-square error of the FDTD results with respect to the Berreman method is ~ 1×10^{-4}.

#### 3.2. Optical activity

For optically active media, the permittivity tensor *$\tilde{\epsilon}$* should be modified to be:

where *G*≃*n*
_{1}Δ*n*_{c}
(Δ*n*_{c}
is the circular birefringence) [20]. We consider a simple planar slab of optically active media with the following material parameters: *n*
_{1,3}=1.5 and Δ*n*_{c}
=0.05. Again, gradient-index AR coatings are applied at the air interfaces. Light traveling through the medium along the *z*-axis experiences a polarization rotation by *π*Δ*n*_{c}*d*/*λ*, where *d* is the thickness of a medium. Since we are most interested in the polarization properties, the Stokes parameters will be examined.

When the incident light polarization is linear and perpendicular to the structure, the normalized Stokes parameters of the transmitted light are given by

where *ξ*_{c}
=Δ*n*_{c}*d*/*λ*. The results obtained from the FDTD simulation show a very good agreement with analytical solutions, as shown in Fig. 4. The mean-square error is ~3.4×10^{-4}.

#### 3.3. Thin phase gratings

We now analyze a binary diffraction grating in isotropic, dielectric media. Its rectangular index profile is illustrated in Fig. 5 and can be represented by

where *n̄* is the average index of refraction and *δn* is the index modulation. The grating parameters can be captured in the permittivity tensor (*ẽ*) by setting *n*
_{1,2,3}(*x*)=*n*(*x*).

A set of analytic expressions for diffraction efficiencies of binary gratings are available[21]:

where the phase angle *φ* is (*πδnd*)/(*λ* cos *θ*_{g}
); *θ*_{g}
is the propagation angle within the grating, which satisfies *n̄* sin *θ*_{g}
=sin *θ*_{inc}
. One can easily see that Eq. (27) is a even function; *η*±*m* are identical. It should be mentioned that these analytical expressions were derived under several assumptions (such as the paraxial limit) [21]. An alternative method to analyze phase gratings is the rigorous coupled-wave analysis (RCWA) [22, 23], which can provide more accurate solutions for most grating problems.

We modeled a binary grating with λ=20*λ*
_{0}, *n̄*=1.5, *δn*=0.5, and *d*=2*λ*_{0}
. Fig. 5(a) and 5(b) show near-field images captured from FDTD simulations at normal and 30° incidence, respectively. Note that this thin phase grating may have many diffracted orders (up to *m*=±19
^{th}
), which are determined by the grating equation given in Eq. (17). We calculated first order efficiency (*η*+1) in each case. Fig. 6(a) shows the first order efficiency as a function of *dnd*/*λ* when *θ*_{inc}
=0. FDTD results show a good agreement with other solutions. Fresnel losses in the diffraction efficiency are observed in both the FDTD and RCWA results, which is not included in the analytical expressions. Fig. 6(b) shows FDTD results when *θ*_{inc}
=30°. As expected, a good agreement between FDTD and RCWA results is found while the analytical expressions fail. The mean-square error is ~4.8×10^{-2} and ~4.2×10^{-2} for *θ*_{inc}
=0 and 30°, respectively. The rectangular shape of a binary grating may be the main reason for the relatively larger numerical errors than a planar structure (such as a slab) with smoothly varying material properties.

#### 3.4. Polarization gratings

Finally, we apply the split-field FDTD method to analyze anisotropic diffraction gratings known as polarization gratings (PGs) [24, 25]. A periodically patterned, uniaxial birefringence in the grating plane, as shown in Fig. 7, modulates the polarization of light instead of its intensity, which leads to unique diffraction characteristics. The optical symmetry axis follows **n**(*x*)=[sin(*πx*/Λ),cos(*πx*/Λ),0].

The two most compelling features of PG diffraction are the potential for 100% diffraction efficiency into a single order by a thin grating and the polarization selectivity of the first diffraction orders. Theoretical and experimental studies on PGs have been reported by several authors [26, 27, 28]. Recently, we demonstrated high quality LC-based PGs for projection displays as a polarization-independent light modulator using photo-alignment techniques [29]. We also reported preliminary numerical analyses of LC-based PGs in Ref. [30], using the same FDTD method described in Section 2.

A set of analytic expressions for diffraction efficiencies of an infinite PG can be obtained from the Jones matrix analysis as follows [29]:

where *ξ*_{l}
=Δ*n*_{l}*d*/*λ* is the normalized retardation, *η*
_{m=0,±1} is the diffraction efficiency of the *m*^{th}
order, and where we have assumed normal incidence and a thin grating (*Q*=2*πλd*/*n*_{o}
Λ^{2} <1). Unlike conventional thin phase gratings, only three diffraction orders (0- and ±1-orders) are present. While the individual first-orders are highly polarization sensitive (to *S*′_{3}), their sum is independent, just as the zeroth-order is independent of the the incident polarization. In addition, diffraction energy is coupled between the 0-order and *±*_{1}
-orders, depending on the normalized retardation *η*_{l}
; the 0-order reaches a minimum (0%) and the ±1-orders show the maximum diffraction (up to 100%) when
${\xi}_{l}=\frac{1}{2}$
.

To study diffraction properties of PGs we model a periodic anisotropic structure by varying *α* in Eq. (5) along the *x*-direction; *η*
_{1,3}=${n}_{\perp}^{2}$, *η*
_{2}=${n}_{\parallel}^{2}$, *α*(*x*)=*πx*/Λ, and *β*=*γ*=0. The permittivity tensor *$\tilde{\epsilon}$* for the spatially varying anisotropy can be written as

The grating parameters are Λ=20*λ*
_{0}, *n*
_{⊥}=1.5, *n*
_{‖}=1.7, Δ*n*_{l}
=*n*
_{‖}-*n*
_{⊥}=0.2, and *d*=5*λ*
_{0}. We vary the input polarization to see the polarization selectivity of the first orders. In addition, grating structures with different ratios of Λ/*λ*
_{0} and varying incident angles *θ*_{inc}
are simulated in order to study the paraxial limit of PGs. Gradient-index AR coatings are applied at air-PG boundaries and following FDTD results are normalized by an AR-coated dielectric slab.

Fig. 8(a) shows diffraction efficiencies (*η*
_{0} and Σ*η*
_{±1}) as a function of normalized retardation (*ξ*_{l}
=Δ*n*_{l}*d*/*λ*) at normal incidence. A sum of *η*
_{±1} can reach ~100% when *η*_{l}
=0.5, 1.5, 2.5, …. We also checked the polarization state of the diffracted orders. The first orders have orthogonal circular polarizations regardless of the input polarization, while the zeroth order has the same polarization as the input. Fig. 8(b) shows the polarization selectivity of the ±1-orders. As expected, the first order diffraction is sensitive to the ellipticity angle *χ* of the incident polarization. The relationship between *χ* and *S*′_{3} is given by *S*′_{3}=sin(2*χ*). We also note that when the input is circularly polarized, only one of the first orders propagates, and is converted to the opposite sense of handedness; a right-hand circular input is converted to the left-hand circular output and vice versa. Fig. 9(a) and 9(b) show near-field images of a PG captured from FDTD simulations when input polarization is linear and circular, respectively. As expected, an excellent agreement between FDTD and analytical results is found. The mean-square error of the FDTD results for both *η*
_{0} and Σ*η*
_{±1} is ~5.8×10^{-4}.

The expressions of Eq. (28) are strictly valid for normal incidence and small diffraction angles, and an analytic description is not available for the case of oblique incidences or when the ratio Λ/*λ* approaches unity. However, our FDTD tool allows for numerical simulation beyond these limits, and is in fact one of the main motives for its development. For example, we modeled a PG with Λ=10*λ* and
$\xi =\frac{1}{2}$
and calculated diffraction efficiencies at different incident angles from -45° to +45°. We normalized each diffracted order by the total transmitted intensity.

Fig. 10(a) shows normalized efficiencies of the three diffracted orders for incident light with right-hand circular polarization, which results in only one propagating mode (*η*
_{-1}≃1) at normal incidence. However, zeroth order efficiency increases as *θ*_{inc}
increases. It can be explained by the fact that the overall effective retardation varies with the incident angle. We also note that the incident angle response of PG diffraction is slightly asymmetric with respect to *θ*_{inc}
=0. We have further studied of the incident angle response of PGs with different ratios of Λ/*λ*(=2.5, 5, 10). Fig. 10(b) shows the sum of *η*
_{±1} as a function of *θ*_{inc}
for three different PGs with linearly polarized light. The degradation in the diffraction efficiency is faster for PGs with a smaller grating pitch, which can reach the paraxial limit at a small incident angle.

## 4. Conclusion

We have developed an efficient FDTD algorithm for wide-band analysis of periodic anisotropic media. PML and periodic boundary conditions are successfully implemented at oblique incidence using the extended split-field update technique. A vectorial near-to-far field transformation is introduced to accurately describe diffraction for high diffraction angles. We have validated our FDTD algorithm by considering the essential optical properties of four stratified and periodic anisotropic structures as well as compared our results to well-known solutions: the twisted-nematic liquid crystal, an optically active slab, the thin binary phase grating, and a polarization grating. Excellent agreement was found in all cases. The FDTD method described in this paper is a useful numerical tool to analyze and design optical elements in periodic (and stratified) media, such as anisotropic gratings and photonic band-gap structures. A package of computer code in standard C/C^{++} format will be made available at http://www.ece.ncsu.edu/oleg/tools as Open Source software, with the name “WOLFSIM (Wideband OpticaL Fdtd SIMulator).”

## Acknowledgements

The authors gratefully acknowledge support from the National Science Foundation (grants ECCS-0621906 and IIP-0539552) and a partnership with ImagineOptix Inc. and Southeast TechInventures Inc.

## References and links

**1. **R. C. Jones, “A new calculus for the treatment of optical systems. I. Description and discussion of the calculus,” J. Opt. Soc. Am. **31**, 488–493 (1941). [CrossRef]

**2. **A. Lien, “Extended Jones matrix presentation for the twisted nematic liquid-crystal display at oblique incidence,” Appl. Phys. Lett. **57**, 2767–2769 (1990). [CrossRef]

**3. **D. W. Berreman, “Optics in stratified and anisotropic media: 4×4-matrix formulation,” J. Opt. Soc. Am. **62**, 502–510 (1972). [CrossRef]

**4. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**5. **E. K. Miller, L. Medgyesi-Mitschang, and E. H. Newman, *Computational electrodynamics - frequency-domain method of moments* (IEEE Press, 1992).

**6. **S. G. Garcia, T. M. HungBao, R. G. Martin, and B. G. Olmedo, “On application of finite methods in time domain to anisotropic dielectric waveguides,” IEEE Trans. Microwave Theory Tech. **44**, 2195–2206 (1996). [CrossRef]

**7. **Y. A. Kao, “Finite-difference time-domian modeling of oblique incidence scattering from periodic surfaces,” Master’s thesis, Massachusetts Institute of Technology, Cambridge, MA (1997).

**8. **A. Taflove and S. C. Hagness, *Computational electrodynamics: finite-difference time-domain method*, 2nd ed. (Artech House, Norwood, MA, 2000), Chap. 13.

**9. **J. A. Roden, S. D. Gedney, M. P. Kesler, J. G. Maloney, and P. H. Harms, “Time-domain analysis of periodic structures at oblique incidence: orthorgonal and nonorthogonal FDTD implementation,” IEEE Trans. Microwave Theory Tech. **46**, 420–427 (1998). [CrossRef]

**10. **S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antennas Propag. **44**, 1630–1639 (1996). [CrossRef]

**11. **J. A. Kong, *Theory of electromagnetic waves* (Wiley, New York, 1975).

**12. **G. B. Arfken and H. J. Weber, *Mathematical methods for physicists*, 4th ed. (Academic Press, San Diego, 1995).

**13. **J. Schneider and S. Hudson, “The finite-difference time-domain method applied to anisotropic material,” IEEE Trans. Antennas Propag. **41**, 994–999 (1993). [CrossRef]

**14. **C. M. Titus, P. J. Bos, J. R. Kelly, and E. C. Gartland, “Comparison of analytical calculations to finite-difference time-domain simulations of one-dimensional spatially varying anisotropic liquid crystal structures,” Jpn. J. Appl. Phys. **38**, 1488–1494 (1999). [CrossRef]

**15. **X. Zhang, J. Fang, K. K. Mei, and Y. Liu, “Calculations of the dispersive characteristics of microstrips by the time-domain finite difference method,” IEEE Trans. Microwave Theory Tech. **36**, 263–267 (1988). [CrossRef]

**16. **C. M. Furse, S. P. Mathur, and O. P. Gandhi, “Improvements to the finite-difference time-domain method for calculating the radar cross section of a perfectly conducting target,” IEEE Trans. Microwave Theory Tech. **38**, 919–927 (1990). [CrossRef]

**17. **P. Yeh and C. Gu, *Optics of liquid crystal displays* (Wiley, New York, 1999).

**18. **C. H. Gooch and H. A. Tarry, “The optical properties of twisted nematic liquid crystal structures with twist angles ≤90 degrees,” J. Phys. D: Appl. Phys. **8**, 1575–1584 (1975). [CrossRef]

**19. **W. H. Southwell, “Gradient-index antireflection coatings,” Opt. Lett. **8**, 584–586 (1983). [CrossRef] [PubMed]

**20. **G. R. Fowles, *Introduction to modern optics* (Holt, Rinehart and Winston, New York, 1968).

**21. **I. Richter, Z. Ryzi, and P. Fiala, “Analysis of binary diffraction gratings: comparison of different approaches,” J. Mod. Opt. **16**, 1915–1917 (1991).

**22. **M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. **71**, 811–818 (1981). [CrossRef]

**23. **M. G. Moharam, E. B. Grann, D. A. Pommet, and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” Proc. IEEE **73**, 894–937 (1985).

**24. **S. D. Kakichashvili, “Method of recording phase polarization holograms,” Sov. J. Quant. Electron. **4**, 795–798 (1974). [CrossRef]

**25. **L. Nikolova and T. Todorov, “Diffraction efficiency and selectivity of polarization holographic recording,” Opt. Acta **31**, 579–588 (1984). [CrossRef]

**26. **I. Naydenova, L. Nikolova, T. Todorov, N. Holme, P. Ramanujam, and S. Hvilsted, “Diffraction from polarization holographic gratings with surface relief in side-chain azobenzene polyesters,” J. Opt. Soc. Am. B **15**, 1257–1265 (1998). [CrossRef]

**27. **J. Tervo and J. Turunen, “Paraxial-domain diffractive elements with 100% efficieincy based on polarization gratings,” Opt. Lett. **25**, 785–786 (2000). [CrossRef]

**28. **G. P. Crawford, J. N. Eakin, M. D. Radcliffe, A. Callan-Jones, and R. A. Pelcovits, “Liquid-crystal diffraction gratings using polarization holography alignment techniques,” J. Appl. Phys. **98**, 123,102 (2005). [CrossRef]

**29. **M. J. Escuti and W. M. Jones, “Polarization independent switching with high contrast from a liquid crystal polarization grating,” SID Digest **37**, 1443–1446 (2006). [CrossRef]

**30. **C. Oh, R. Komanduri, and M. J. Escuti, “FDTD and elastic continuum analysis of the liquid crystal polarization grating,” SID Digest **37**, 844–847 (2006). [CrossRef]