5. Diffraction; Transforming FDTD Data into Far Field, for Modeling Microdisplays and Gratings

6.  Some Results.
 

1. Introduction:  One-D vs. Two-D Computer Simulation Methods (Chuck Titus)
 

In the past, optical simulation of liquid crystal structures has been dominated by “one-dimensional” methods, such as the extended Jones method[A] and the Berreman 4 X 4 matrix method[B].  These methods are restricted to or designed specifically for simulations of inhomogeneous liquid crystal (anisotropic) structures with variations of the liquid crystal director in only one direction, perpendicular to the glass substrates enclosing the liquid crystal.
 

All pixel areas in liquid crystal displays include at least some regions in which the liquid crystal director (i.e. the dielectric tensor) also varies in a direction parallel to the glass substrates.  Such variations are caused by effects such as disclinations (e.g. reverse tilt in a TN) and fringe field areas near the edges of each pixel (i.e. near the edges of the electrodes which define each pixel).  However, in a low-resolution direct-view display, the pixels are large compared to areas covered by disclinations and fringing fields, and thus the contribution of such regions to the total optical performance of a pixel can be ignored.  As a result, the above-mentioned one-dimensional optical simulation methods can be satisfactorily employed for optical analysis of relatively simple low-resolution liquid crystal displays.
 

With the recent emergence of liquid crystal microdisplay devices such as those in projection and head-mounted displays , pixel size is decreasing.   In addition, liquid crystal display structures such as in-plane switching may possess intentionally significant multidimensional spatial inhomogeneities of the liquid crystal director orientation.  Because of this, disclinations and fringing field regions contribute more significantly to net optical performance of each liquid crystal display pixel.  Furthermore, the smaller size of these microdisplay pixels means that diffractive effects also play an increasingly significant role in optical performance.
 

There have been some attempts to use the one-dimensional methods in these cases.  One such technique involves dividing the structure into layers and columns.  Each column is regarded as an independent one-dimensional structure.  One of the one-dimensional methods is then applied to each column.  However, such adaptation of one-dimensional techniques to two-dimensional problems cannot properly take into account lateral scattering of light propagating within the liquid crystal structure.  Nor can those methods accurately calculate the propagation of obliquely incident light through a liquid crystal structure with director variations in more than one dimension.
 

There are other interesting (non-display) liquid crystal devices for which past optical simulations required assumptions and simplifications which impacted the accuracy of the simulation.  One such application involves the use of liquid crystal phase gratings for beam steering applications[C,D].  In those devices, the phase profile is achieved by applying a distribution of voltages to an array of separated electrodes on the surface of a liquid crystal film.  Optical simulations of such structures have been conducted which assume simplified liquid crystal director configurations.  In those analyses, it is assumed that the liquid crystal director varies only in the direction parallel to the grating surface, ignoring the effect of fringing fields.
 

For these and other reasons, optical modeling of such structures may require more robust computational methods than the existing one-dimensional optical simulation techniques.  Among the two-dimensional computational methods, the Finite-Difference Time-Domain (FDTD) method is an approach we are considering.  FDTD makes use of direct discretization of Maxwell’s equations.  Because of this, the FDTD method can accommodate multidimensional inhomogeneity of the dielectric tensor.
 

A)  A. Lien; Appl. Phys. Lett., 57 2767 (1990).

2. How the Finite-Difference Time-Domain (FDTD) Method Works  (Chuck Titus)
 

The FDTD method has been in use in the electrical engineering community for years, but has only recently been applied to liquid crystal structures[E].  FDTD, as introduced by Yee[F] for isotropic media, makes use of direct discretization of the general form of Maxwell’s equations.  Because of this, the FDTD method can accommodate multidimensional inhomogeneity of the dielectric tensor.
 

At the heart of the FDTD method is the spatial discretization of the structure to be analyzed, and the use of central finite differences to approximate Maxwell’s equations.  An oscillating wave source is introduced into the grid and “directed” at the structure to be analyzed.  The entire grid is then time-stepped until steady state response oscillations are obtained at all points on the grid.
 

The previous FDTD simulation of liquid crystals employs this technique, modified for sourceless inhomogeneous anisotropic media[E]:
 

E)  B. Witzigmann, P. Regli, W. Fichtner; J. Opt. Soc. Am. A, 15, 753 (1998).

2.A.  Discretization of Maxwell’s Equations  (Chuck Titus)

Yee[F] introduced the FDTD method for isotropic media.  At the heart of the FDTD method is the discretization of Maxwell’s equations, which can be written, for sourceless inhomogeneous anisotropic media:

(1)

(2)

where e^-1(r) is the inverse of the spatially varying dielectric tensor.  After discretization in time employing the central finite difference, eqs. (1) and (2) can be used as follows to update the three vector components of both fields[E]:
 

(3)

(4)

where subscripts i,j,k,l=(x,y,z) indicate rectilinear coordinate axes, and summation convention is used.  The symbol Î_ijk represents the Levi-Civita permutation symbol.  The i,j component of the inverted dielectric tensor at any grid location r is denoted e^-1_ij(r).  The superscript n denotes the discrete time coordinate, and D t is the time step resulting from finite difference approximation of the time derivative.  Spatial derivatives above are not discretized for the sake of simplicity, but they, too, are approximated with central finite differences.
 

For example, let us look at the case of a linear isotropic dielectric.  In a one-dimensional problem whose grid is the z-axis, the partial derivatives with respect to x and y do not exist.  In that case, full discretization of the x-component of eq. (4) would yield an “update equation for H at the z=(K+½)D h node:
 

(5)

where D h is the separation of the two grid nodes surrounding the location at which H is being updated.  Both the space and time finite difference approximations shown above are second-order accurate.  A graphical illustration of the one-dimensional space-time grid utilized above is shown below.

The entire structure to be analyzed is spatially discretized and its properties are distributed onto the computational grid.  In this grid, electric and magnetic field quantities occupy alternating spatial and temporal locations.  This permits approximation of both space and time derivatives with central finite difference formulae.  This spatial and temporal interleaving, as employed above to update the desired quantities at each grid location, is known as a staggered-leapfrog scheme.
 

E)  B. Witzigmann, P. Regli, W. Fichtner; J. Opt. Soc. Am. A, 15, 753 (1998).

2.B  Illumination of Objects on the FDTD Grid  (Chuck Titus)

Plane Wave.  Illumination of liquid crystal structures is provided by means of the separate-field formulation[G].  In this technique, the inner computational domain is divided into two regions.  In the inner region the total field, the superposition of the source and the scattered fields is calculated.  In the outer region, only the scattered field is calculated.  The oscillating wave source is introduced at the interface between these two regions, added into Maxwell’s equations on the inner side of the interface and subtracted on the outer side.  This technique can produce diffraction-free plane waves.  The separate-field formulation also permits easy extension of this method to off-axis incident light.

Is a plane wave the appropriate type of source?  One drawback of the above plane wave method is diffraction resulting from illumination of a finite object with a plane wave of infinite extent.  Because the finite object has “edges”, the plane wave will diffract about those edges.  In a real liquid crystal display or device, the pixel being simulated is part of a larger structure, so the pixel does not possess the discontinuities suggested in the above figure.
 

We could simply append the appropriate neighboring liquid crystal structures to the left and right of the structure in the above figure.  The virtual source would then also be extended, keeping it outside the new larger structure.  As a result, diffraction resulting from artificial termination of the liquid crystal structure on the computational grid would be pushed away from the region of interest.
 

One previous FDTD simulation of liquid crystals[E] employed periodic boundaries (left and right sides in the figure above), in effect assuming an infinite lateral juxtaposition of identical liquid crystal structures.  In this case the structure is also illuminated by what is, in effect, a plane wave of infinite extent.  This may not always be desirable.
 

There is something not quite proper about the use of a plane wave.  For instance, one may wish to simulate the performance of a liquid crystal grating.  Such a grating may be used for steering an incident laser beam into a particular direction.  The laser is a “plane wave” only at points where a minimum beam waist is achieved.  In general, a propagating laser beam is not a plane wave, and its intensity is not spatially uniform.
 

The near field results obtained by probing a liquid crystal structure with a plane wave are not, in general, identical to those obtained by probing the same structure with a TEM₀₀ laser beam.  In addition, computed far field results from an array of pixels illuminated by a plane wave of uniform intensity will differ from those obtained when illuminating the same structure with a gaussian beam.  One can “fix” this by applying a gaussian weighting to the transmitted fields resulting from plane wave illumination, but that still does not accurately model the true probe-and-response.
 

In order to eliminate the above concerns, a laser source has been implemented.  More details will be forthcoming.
 

E)  B. Witzigmann, P. Regli, W. Fichtner; J. Opt. Soc. Am. A, 15, 753 (1998).

2.C. Grid Termination, Boundary Conditions  (Chuck Titus)
 

In order to produce an accurate solution to Maxwell’s equations, the results must be a solution over an infinite space.  However, practical limitations of computer power and memory require the termination of the computational grid.  Any such termination method must not affect the computations inside the finite computational grid.  In other words, the grid termination must not result in reflection of propagating waves back into the grid.
 
 

There are a variety of grid termination techniques[G], but only two are discussed here.  The early FDTD liquid crystal study employed the second-order Mur treatment[H]for terminating the grid.  The Mur grid termination technique applies a numerical solution to the wave equation along the grid boundary.  It is computationally inexpensive and relatively easy to implement.
 
 

Our testing of the Mur absorbing boundary revealed that it did not result in efficient termination except for waves normally incident (or nearly so) upon the grid boundary.  Reflections of waves incident upon the Mur boundary were, over many time steps, of sufficient magnitude to influence FDTD results.  For this reason, the more effective “Perfectly matched layer” (PML) absorbing layer[I] is often more useful.
 

The PML is implemented as an absorbing layer which surrounds the inner computational domain.  It consists of a fictitious lossy medium which incorporates magnetic as well as electric conductivity:
 

(6)

(7)

Impedance matching of the PML with the adjacent inner computational domain makes that interface very transmissive to waves radiating from the inner grid.  As a result, radiation reaching the boundary passes into the PML without reflection.  The electric and magnetic conductivities of the PML then result in rapid decay of waves passing into it.
 

G)  A. Taflove; Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston (1995).

H)  G. Mur; IEEE Trans. Electr. Compat. 23, 377, (1981).

I)   J-P. Berenger; J. Comp. Phys., 114, 185 (1994).

2.D. Numerical Stability and Accuracy  (Chuck Titus)
 

Stability:  Two other matters are of importance in FDTD simulations.  First, limiting the size of the time step is necessary to prevent growth of numerically induced oscillations.  Second, inaccuracies can result from the temporal and spatial discretization of the problem.
 

Limiting the size of the time step is necessary in order to prevent growth of numerically induced oscillations.  The limit for an N-dimensional square FDTD grid employing second-order accurate spatial and temporal finite differences is expressed in the following relation[J]:
 

(8)

where D t is the size of the time step and v_max is the maximum phase velocity of waves traveling on the grid.  The spatial step  D  h is the size of the standard Yee cell (two node separations across), normally defined as  l  /N, where N is some integer.
 
 

Accuracy:  FDTD inaccuracies result from the temporal and spatial discretization of the problem.  For instance, the second-order accurate central finite difference approximations used in this study are somewhat less than exact.  Such discretization can manifest itself in numerical dispersion [G] , which is the departure of numerical phase velocity from that which would occur in continuous space.  In particular, the phase velocity of waves in an FDTD simulation employing second-order accurate spatial and temporal finite differences is always less than the continuous-space value.  In other words, the numerical wavelength is always less than the continuous-space value, and spectral data generated by FDTD simulations will be shifted to longer wavelengths
 

The magnitude of this shift depends on spacing between nodes on the grid, time step size, dielectric properties, and the frequency of the propagating wave[K].  Because time-stepping to a steady-state FDTD result may require accounting for multiple reflections, the impact of numerical dispersion also depends on the size of the structure to be analyzed and the number of time steps required to reach steady state.
 

Of these factors, only the time step size and node spacing are not predetermined by the problem at hand.  One might expect that as time step and node spacing approach zero, errors are reduced and the FDTD result approaches the exact solution.  This behavior for node spacing has been confirmed[K].  However, as has been found in previous studies[L,M], FDTD results are more accurate when the time step is equal to, not less than, its maximum stable value.  This may be due in part to the fact that as the time step is decreased, more steps are necessary to reach steady state, increasing the impact of roundoff error.
 

In addition, the continuous spatial variations of the dielectric tensor of liquid crystal structures are discretized.  This may lead to field discontinuities and spurious reflections which are numerical in origin and not at all representative of actual light propagation through such structures.  We have not yet determined if this is a significant source of error.  If such error does occur, it will decrease as the node spacing is decreased.
 

Thus, for a given discretization scheme, the only simple control over accuracy at our disposal is the choice of node spacing.
 

G)  A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston (1995).

J)  A. Taflove, M. E. Brodwin:  IEEE Trans. Microwave Theory &Techniques 23, 623, (1975).

K)  C. M. Titus, P. J. Bos, J. R. Kelly, E. C. Gartland, Japanese Journal of Applied Physics, Part 1, 38, 1488, (1999).

L)  I. S. Kim, W. J. R. Hoefer:  1989 IEEE AP-S Symp. Dig., 1108, (1989).

M)  I. S. Kim, W. J. R. Hoefer:  Electron. Lett. 26, 485, (1990).

2.E. Computational Expense  (Chuck Titus)
 

Of course, as node spacing is reduced for a fixed geometry FDTD problem, memory requirements increase and more time steps are required to reach steady state.  For instance, if the node spacing is cut in half, a two-dimensional problem must be represented by four times as many nodes.  And since the maximum stable time step size is proportional to node spacing, cutting the node spacing in half also necessitates twice as many time steps to reach a steady state solution.  In general, computational hardware expense increases roughly as the square of the relative increase of number of nodes per wavelength, and computational time increases roughly as the cube of the relative increase of number of nodes per wavelength.
 

The following is an illustration the hardware required for a specific liquid crystal problem.

As will be shown below, use of fourth-order accurate finite differences[G] reduces the required grid spacing to one-twentieth of a wavelength, or one-ninth the number of nodes on a two-dimensional grid.  The increased computational load brought on by use of fourth-order finite differences is more than offset by the accompanying reduction in the number of grid points.
 

Additional time savings can be obtained by parallellizing the computational code and executing it on a multi-processor system.  This has just been completed.  More details later.
 

G)  A. Taflove; Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston (1995).

3. FDTD Capabilities (Chuck Titus)
 

The FDTD method provides an avenue toward optical simulation of liquid crystal devices with small pixel size and/or multi-dimensional inhomogeneities of director orientation.  These conditions occur in liquid crystal microdisplays, some diffractive liquid crystal devices, and in existing liquid crystal displays such as those employing in-plane switching.  In future comparisons with far-field results of one-dimensional optical techniques applied to these problems, the one-dimensional methods may well prove equal to the task.  However, without very precise experimental methods or a tool such as FDTD, that cannot be known.
 

The FDTD method, as described in previous sections, is capable of directing laser beams or plane waves at arbitrary incidence onto any liquid crystal structure placed on the grid.  After numerical time-stepping to steady state, transmitted and reflected responses can be collected.  Those responses can be utilized to compute transmitted and reflected intensities along the surface of a liquid crystal pixel or transformed into the far field (i.e. diffraction).
 

As a further illustration of the ability of FDTD to simulate a variety of optical problems, the figure below shows the interference pattern resulting from two oscillating dipoles placed on an FDTD grid.

4. FDTD vs. Analytical Solutions to Liquid Crystal Problems  (Chuck Titus)
 

It must not be forgotten that, as with any discretization of a continuous physical model, the FDTD method provides only a numerical approximation of the solution to Maxwell’s equations.  The accuracy of this method has been tested for a variety of structures and media, including some anisotropic media[N].  None of these tests, however, cover the inhomogeneous dielectric anisotropy characteristic of advanced liquid crystal displays.  In particular, liquid crystal structures can contain disclinations which result in large variations of the dielectric tensor over sub-wavelength distances.  Some analysis is necessary in order to know how to apply FDTD modeling to liquid crystal structures in cases ranging from slowly varying director orientation found in in-plane switching to structures containing defects.
 

Toward that end, we reported a comparison of FDTD and analytic solutions for several well-known liquid crystal problems with exact solutions[K].  One of the test cases was Bragg reflection from a planar cholesteric layer.  Uniform cholesteric liquid crystal layers are known to exhibit Bragg reflection of wavelengths comparable to the cholesteric pitch.  An exact analytical solution for this problem at normal incidence can be found in reference P.
 

As spatial variation of the dielectric tensor increases, finer node spacing is required in order to obtain accurate FDTD computation.  It is not the spatial rate of change of the liquid crystal director that is important.  Maxwell’s equations do not explicitly consider the orientation of molecules.  Rather, it is the spatial rate of change of the dielectric tensor, which is influenced by both the molecular orientation and the molecular dielectric anisotropy.
 

How fine the grid should be will depend on the significance of defects in the liquid crystal structure being modeled.
 

All derivatives in our report[K] were approximated with second-order-accurate finite differences.  The study concluded that, without resorting to higher-order finite-differences, choice of node spacing was the only simple means of achieving accurate FDTD results for liquid crystal structures.  Numerical dispersion caused significant shift of numerical results to longer wavelengths unless an extremely fine grid was employed.  The use of such a fine grid in a two-dimensional computation would not have been practical, requiring more memory than could be installed on most existing desktop computers.  In addition, very long computation times would have been required.
 

Since that study, this problem has been addressed by replacing the second-order finite difference approximation of the spatial derivatives with a fourth-order approximation.  FDTD computation of the same cholesteric layer produced accurate results for l /20 node spacing (see figure below).

K)  C. M. Titus, P. J. Bos, J. R. Kelly, E. C. Gartland; Japanese Journal of Applied Physics, Part 1, 38, 1488, (1999).

N)  J. Schneider, S. Hudson; IEEE Trans. Antenna. & Propag. 14, 994, (1993) .

P)  H. Wohler, G. Haas, M. Fritsch, D. A. Mlynski; J. Opt. Soc. Am. A, 5, 1554 (1988).
 
 
 

The FDTD method described above is used to direct a light source at a liquid crystal structure and then compute the electromagnetic fields exiting that structure.  Those exit fields are in the near field region, which is defined as being within one wavelength or so beyond the exit surface.  The next step in the simulation is to utilize those results to obtain the distribution of light far from the liquid crystal structure.  That far zone distribution of light is what we often call the diffraction pattern.  In the case of a beam steering liquid crystal grating, the diffraction pattern includes the portion of incident light which is steered into the desired direction, plus that portion which propagates in other directions.  Complete knowledge of this distribution of diffracted light may be of use in improving the design of these liquid crystal beam steering devices.
 

After employing the FDTD method to obtain the near fields, the same technique could be applied outward through space to some distance from the diffracting liquid crystal structure.  Assuming careful attention is paid to numerical issues, one could thusly obtain the far field pattern.  However, as has already been discussed, the FDTD method incurs substantial computational costs even for relatively small structures.  Because of limits to computer memory and speed, direct computation of far field behavior is not feasible at this time.  An alternate technique is required.  For this we turn to diffraction theory.
 

The general topic of diffraction theory is discussed in more detail in another section.  The following is the paraphrased conclusion of that document, which is derived from a detailed study of diffraction theory.
 

The Fraunhofer approximation is just that, an approximation.  First of all, its derivation assumes the distance from diffracting object to observing screen is large compared to the size of the diffracting object and a wavelength of light.  Secondly, and not as commonly known, the Fraunhofer approximation produces results which are increasingly inaccurate as the diffraction angle increases.  This second concern has led us to the use of the more accurate Kirchhoff surface integral.
 

This has also been implemented in code, and is used to compute the far field diffraction pattern resulting from laser illumination of two-dimensionally inhomogeneous liquid crystal structures.  The accuracy of FDTD near field computation plus Kirchhof surface integral far field has been verified by comparing computed far field propagation of gaussian beams with theoretically-predicted values.  Other checks were performed, computing refraction of gaussian beams by isotropic prisms and comparing with Snell’s law.  Additional checks included comparison of computed diffraction of gaussian beams by ideal blazed gratings with theory.  Those results will appear here in the future.

6. Some Results. (Chuck Titus)
 

The following data are results produced by the use of FDTD to compute near fields transmitted through liquid crystal structures, and subsequent transformation into the far field by means of the Kirchhoff surface integral.  These are preliminary results, prior to implementation of the laser source.
 

The first structure is a twisted-nematic (TN) pixel.  Voltage above the Freedericksz Threshold, applied to the TN structure, results in a liquid crystal structure which includes a reverse-tilt disclination. The two-dimensional liquid crystal director map, obtained courtesy Dr. Jack Kelly, appears as:

This structure was probed with a plane wave, incident from below.  The plane wave source was linearly polarized perpendicular to the plane of the structure.  FDTD computation yielded the surface intensity distribution shown below (left side).  The far field computation produced the second plot (to the right).