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.

 

Fig. 1. Schematic views of (a) a 2D FDTD computational space and (b) a bi-dimensional FDTD grid, where field variables locate.

Download Full Size | PDF

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, ε ₀ and µ ₀ 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 ₁,e ₂,e ₃) 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

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² θ_inc , B=κ_yz sin³ θ _inc , C=(1/A)κ_yz sin² θ_inc , and D=1 - (B/A)κ_zy sin θ_inc -κ_yy sin² θ_inc . The discrete version of Eq. (12a) can be found:

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 ₀=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 ₀, 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 ω ₀ is the peak-frequency, t ₀ 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 λ ₀ is either the wavelength for monochromatic light or the peak-wavelength of a Gaussian pulse corresponding to ω ₀.

Table 1. FDTD simulation parameters

View Table

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_ld/λ, 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_ld/λ. When the twist angle is 90°, the maximum transmittance occurs when:

 

Fig. 2. A 90° twisted-nematic liquid crystal cell between cross-polarizers. Red arrows depict the nematic director of liquid crystals and d is the cell thickness.

Download Full Size | PDF

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$, ε ₂=$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 ₁ is the refractive index of the incoming media, n ₂ is the index for the outgoing media, and t is the AR thickness; we set t=λ ₀ 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 ₀ where S ₀ 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. 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.

 

Fig. 3. FDTD results of a 90° twisted-nematic liquid crystal cell between ideal cross-polarizers: (a) the polarization state right after the LC layer; (b) the transmittance (T) through the TN-LC cell with cross polarizers. LC material parameters are n _⊥=1.5 and Δn_l =0.2. Gradient-index AR coatings are applied at the air-LC boundaries and the inset figure shows the transmittance of an AR-coated dielectric slab with n=1.6.

Download Full Size | PDF

3.2. Optical activity

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

where G≃n ₁Δ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_cd/λ, where d is the thickness of a medium. Since we are most interested in the polarization properties, the Stokes parameters will be examined.

 

Fig. 4. FDTD results of the normalized Stokes parameters for the transmitted light from a slab of optically active media with circular birefringence Δn_c =0.05.

Download Full Size | PDF

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_cd/λ. 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]:

 

Fig. 5. Near-field images of a binary phase grating for (a) θ_inc =0 and (b) θ_inc =30°. Grating parameters are Λ=20λ, n̄=1.5, δn=0.5, and d=2λ.

Download Full Size | PDF

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λ ₀, n̄=1.5, δn=0.5, and d=2λ₀ . 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.

 

Fig. 6. FDTD results of the first order efficiency (η+1) of a binary phase grating at (a) θ_inc =0 and (b) θ_inc =30°. Grating parameters are Λ=20λ_o , n̄=1.5, δn=0.5, and d=2λ_o .

Download Full Size | PDF

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.

 

Fig. 7. Periodically varying anisotropy profile of a polarization grating; red arrows depict induced polarization state by local anisotropy. Λ is the grating pitch.

Download Full Size | PDF

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_ld/λ 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 Λ² <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′₃), 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 ±₁ -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$, η ₂=$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λ ₀, n _⊥=1.5, n _‖=1.7, Δn_l =n _‖-n _⊥=0.2, and d=5λ ₀. We vary the input polarization to see the polarization selectivity of the first orders. In addition, grating structures with different ratios of Λ/λ ₀ 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 (η ₀ and Ση _±1) as a function of normalized retardation (ξ_l =Δn_ld/λ) 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′₃ is given by S′₃=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 η ₀ and Ση _±1 is ~5.8×10^-4.

 

Fig. 8. FDTD results of diffraction efficiencies of a linear PG as a function of (a) normalized retardation Δn_ld/λ and (b) ellipticity angle (χ) of the polarization ellipse. Grating parameters are Λ=20λ ₀, d=5λ ₀, n _⊥=1.5, and Δn_l =0.2.

Download Full Size | PDF

 

Fig. 9. Near-field images of the PG for different input polarizations: (a) arbitrary linear input polarization; (b) circular input polarization. Grating parameters are Λ=20λ, n _⊥=1.5, Δn_l =0.2 and $\frac{\Delta n_{l}d}{\lambda }=\frac{1}{2}$ .

Download Full Size | PDF

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.

 

Fig. 10. FDTD results of the incident angle response of PG diffraction for: (a) right-hand circularly polarized light; (b) three different grating pitches (Λ/λ=2.5,5, 10). Grating parameters are n _⊥=1.5, Δn_l =0.2, and d=2.5λ.

Download Full Size | PDF

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).”