OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 17, Iss. 13 — Jun. 22, 2009
  • pp: 10895–10909

Optics InfoBase > Optics Express > Volume 17 > Issue 13 > Page 10895

« Show journal navigation

A finite element beam propagation method for simulation of liquid crystal devices

Pieter J. M. Vanbrabant, Jeroen Beeckman, Kristiaan Neyts, Richard James, and F. Anibal Fernandez  »View Author Affiliations


Optics Express, Vol. 17, Issue 13, pp. 10895-10909 (2009)
http://dx.doi.org/10.1364/OE.17.010895


View Full Text Article

Acrobat PDF (2928 KB)


Advanced Search

Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

An efficient full-vectorial finite element beam propagation method is presented that uses higher order vector elements to calculate the wide angle propagation of an optical field through inhomogeneous, anisotropic optical materials such as liquid crystals. The full dielectric permittivity tensor is considered in solving Maxwell’s equations. The wide applicability of the method is illustrated with different examples: the propagation of a laser beam in a uniaxial medium, the tunability of a directional coupler based on liquid crystals and the near-field diffraction of a plane wave in a structure containing micrometer scale variations in the transverse refractive index, similar to the pixels of a spatial light modulator.

© 2009 Optical Society of America

1. Introduction

The beam propagation method (BPM) is a versatile and efficient numerical method, well suited to the analysis of light propagation in today’s complex photonic devices. This is reflected in the large variety of BPM algorithms that have been studied, ranging from methods based on fast fourier transforms (FFT) [1

J. M. López-Doña, J. G. Wangüemert-Pérez, and I. Molina-Fernández, “Fast-fourier-based three-dimensional full-vectorial beam propagation method,” IEEE Photonics Technol. Lett. 17, 2319–2321 (2005). [CrossRef]

], finite difference (FD) schemes [2

C. Ma and E. Van Keuren, “A three-dimensional wide-angle BPM for optical waveguide structures,” Opt. Express 15, 402–407 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?uri=OE-15-2-402. [CrossRef] [PubMed]

, 3

A. J. Davidson and S. J. Elston, “Three-dimensional beam propagation model for the optical path of light through a nematic liquid crystal,” J. Mod. Opt. 53, 979–989 (2006). [CrossRef]

, 4

Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. 23, 2014–2019 (2006). [CrossRef]

] and finite element (FE) [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

, 6

D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, “Mixed finite element beam propagation method,” J. Lightwave Technol. 16, 1336–1342 (1998). [CrossRef]

] implementations. Employing a finite element scheme offers high accuracy and the flexibility to model arbitrarily curved structures. Hybrid edge/nodal finite elements [7

M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). [CrossRef]

] have been used with great success in modeling vector electric or magnetic fields. These elements exclude spurious solutions whilst incorporating the desired continuity conditions at dielectric interfaces.

Our goal is the development of a full-vector finite element BPM to simulate the evolution of the electric field of an input optical field upon propagation through inhomogeneous anisotropic media. In particular, we want to apply the method for optical modeling of liquid crystal devices such as liquid crystal displays (LCDs), spatial light modulators (SLMs) or tunable photonic components based on liquid crystals. In order to represent the inhomogeneous orientation of the anisotropic liquid crystal, the full dielectric tensor must be considered in the beam propagation method, as described in Section 2. Furthermore, a wide angle propagation scheme is desired without first order approximation of the second order derivative that appears in the wave equation for the electric field. Perfectly matched layers (PMLs) for anisotropic materials [8

F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microw. Guid. Wave Lett. 8, 223–225 (1998). [CrossRef]

, 5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

] are positioned at the edges of the modeling window to enforce absorbing boundary conditions in order to avoid nonphysical reflection of the considered optical field at the edges.

To demonstrate the wide applicability of the presented beam propagation algorithm, we illustrate in Section 3 how the method can be used to simulate the light propagation in three applications involving anisotropic materials of arbitrary orientation. A first example illustrates, as a test example, the propagation of a linearly polarized gaussian light beam through a homogeneous uniaxial material. Next, the tunability of the coupling characteristics of a directional coupler achieved by changing the orientation of a liquid crystal cladding layer is illustrated to show the long distance stability of the developed BPM. In the third example the diffraction of a plane wave propagating through a structure containing micrometer variations in permittivity, such as the pixels in a spatial light modulator (SLM), is discussed and the influence of the pixel dimensions on the near-field light diffraction is shown.

2. Theory: the beam propagation method for materials with general anisotropy

2.1. Finite element discretization of the wave equation

The starting point to consider the evolution of time harmonic fields is the Helmholtz equation for the electric field Ē(x,y, z), as derived from the Maxwell equations:

∇× ( μ̿ −1∇× Eˉ) − k02 ε̿. Eˉ=0,
(1)

where k 0 is the wavenumber in vacuum, ε̿ is the relative permittivity tensor and µ̿ is the relative permeability tensor. Perfectly matched layers (PML) for anisotropicmedia [8

F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microw. Guid. Wave Lett. 8, 223–225 (1998). [CrossRef]

] for reflectionless absorption of electromagnetic waves can be incorporated in Eq. (1) by interchanging µ̿−1 and ε̿ with p̿ and q̿ respectively. These are permittivities and permeabilities modified to include the absorption coefficients in the eight PML regions [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

] at the borders of the computationalwindow, as shown in Fig. 1(a). The directional coupler shown in Fig. 1(a) will be further considered for BPM analysis in Subsection 3.2.

Fig. 1. (a) Transverse cross section of a directional coupler, indicating the different PML regions 1–8 at the edges of the computational window, (b) Sketch of a hybrid linear tangential/ quadratic normal (LT/QN) vector element indicating the transverse and nodal unknowns Φ ti (i=1 to 8) and Φ zj (j=1 to 6), respectively.

Using a slowly varying envelope approximation, the electric field Ē can be separated into a slowly varying complex field Φ¯ (x,y, z)=ϕx (x,y, z)1̄ x +ϕy (x,y, z)1̄ y +ϕz (x,y, z)1̄ z and a phase factor exp(−jk 0 n 0 z), where z is the direction of light propagation:

Eˉ ( x,y,z)= Φˉ ( x,y,z)exp ( −j k0 n0z),
(2)

with n 0 an appropriate reference refractive index. In the finite element method, Eq. (1) is discretized by dividing the structure cross section Ω into small hybrid edge/nodal triangular finite elements. Within each element, the vector transverse and longitudinal field components are expanded as:

[ ϕx ϕy ϕz]= [ [ Nx]T [0] [ Ny]T [0] [0] j[L]T] [ ϕte ϕze],
(3)

where ϕe t and ϕe z are the edge and nodal values, respectively, in the element being considered. Hybrid LT/QN elements are chosen, as shown in Fig. 1(b), that use linear tangential (LT) Nx and Ny and quadratic normal (QN) L shape functions for interpolation of the transverse and longitudinal field [7

M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). [CrossRef]

]. The use of these vector elements avoids the occurrence of spurious solutions and leads to a simple and correct description of the continuity condition of the electric field at dielectric interfaces. LT/QN shape functions are chosen over first order functions due to their higher-order convergence. The full dielectric permittivity tensor ε̿ is maintained and a Galerkin procedure [9

J. Jin, The finite element method in electromagnetics , 2nd edition (Wiley, New York US, 2002).

] is applied to the wave equation from Eq. (1). Substituting Eqs. (2) and (3) into Eq. (1) leads to the beam propagation expression:

[ [ Btt] [0] [0] [0]] d2 dz2 { [ ϕt ϕz]}− [ 2j k0 n0 [ Btt] j[ Btz] j [ Bzt] [0]] ddz { [ ϕt ϕz]}

− { [ [ Att] −j[ Ctz] j[ Czt] [ Bzz]] −j k0 n0 [ [0] j[ Btz] j [ Bzt] [0]]+ k02 n02 [ [ Btt] [0] [0] [0]]} { [ [ ϕt] [ ϕz]]}=0
(4)

For structures invariant with respect to the propagation direction, this equation may be arranged as an eigenvalue problem for the modes of fully anisotropic liquid crystal waveguides [10

J. Beeckman, R. James, F. A. Fernandez, W. De Cort, P. J. M. Vanbrabant, and K. Neyts, “Calculation of fully anisotropic liquid crystal waveguide modes,” accepted for publication in J. Lightwave Technol. (2009). [CrossRef]

]. Assuming the variation in ∂pii/∂z is sufficiently small, the submatrices in Eq. (4) can be evaluated as:

[ Att]= Σe ∬e pzz ( ∂{ Nx} ∂y ∂{ Nx}T ∂y+ ∂{ Ny} ∂x ∂{ Ny}T ∂x − ∂{ Nx} ∂y ∂ { Ny}T ∂x − ∂{ Ny} ∂x ∂ { Nx}T ∂y)

− k02 ( qxx { Nx} { Nx}T + qxy { Nx} { Ny}T + qyx { Ny} { Nx}T + qyy { Ny} { Ny}T)dxdy
(5)
[ Btt]= Σe ∬e [ pyy { Nx} { Nx}T+ pxx{ Ny} { Ny}T]dxdy
(6)
[ Btz]= Σe ∬e [ pyy { Nx} ∂ {L}T ∂x+ pxx { Ny} ∂ {L}T ∂y]dxdy
(7)
[ Btz]= Σe ∬e [ pyy ∂ {L} ∂x { Nx}T+ pxx ∂ {L} ∂y { Ny}T]dxdy
(8)
[ Bzz]= Σe ∬e [ pyy ∂ {L} ∂x ∂ {L}T ∂x+ pxx ∂ {L} ∂y ∂ {L}T ∂y − k02 qzz {L} {L}T]dxdy
(9)
[ Ctz]= Σe ∬e k02 [ qxz { Nx} {L}T+ qyz{ Ny} {L}T]dxdy
(10)
[ Czt]= Σe ∬e k02 [ qzx {L} { Nx}T+ qzy {L} { Ny}T]dxdy,
(11)

where the coefficients related to the transverse anisotropy in permittivity (qxy and qyx ) appear in [Att ] and the coefficients related to the longitudinal anisotropy in permittivity (qxz, qzx, qyz and qzy ) appear in [Ctz ] and [Czt ]. The ϕz component can be eliminated from the second row of Eq. (4):

ϕz=−j Bzz −1 Bzt d ϕtdz − Bzz −1 ( j Czt+ k0 n0 Bzt) ϕt.
(12)

Substituting the expression for ϕz into the first row of Eq. (4) yields a second order differential equation for the transverse field component:

{ [ Btt]− [ Btz] [ Bzz] −1 [ Bzt]} d2 ϕt dz2 − {2j k0 n0 [ Btt] −j [ Btz] [ Bzz] −1 (j [ Czt] + k0 n0 [ Bzt]) +j (j [ Ctz] − k0 n0 [ Bzz]) [ Bzz] −1 [ Bzt] dϕtdz} − { [ Att]+ k02 n02 [ Btt] + (j [ Ctz] − k0 n0 [ Btz]) [ Bzz] −1 (j [ Czt]+ k0 n0 [ Bzt])} ϕt=0,
(13)

which is the equation that governs the wide angle beam propagation of the transverse field.

2.2. Wide angle beam propagation algorithm

In order to derive a recurrence scheme for the light propagation in the z-direction, Eq. (13) is formally rewritten in the form of a first order differential equation:

A ddz [ ϕt dϕtdz]+B [ ϕt dϕtdz] =0,
(14)

where A and B are defined as:

A= [ [ A11] [ A12] [1] [0]]
(15)
B= [ [ B11] [0] [0] − [1]],
(16)

with [1

J. M. López-Doña, J. G. Wangüemert-Pérez, and I. Molina-Fernández, “Fast-fourier-based three-dimensional full-vectorial beam propagation method,” IEEE Photonics Technol. Lett. 17, 2319–2321 (2005). [CrossRef]

] the identity matrix having the same dimensions as the submatrices [A 11], [A 12] and [B 11] which are defined as:

[ A11]= −2j k0 n0 [ Btt] − [ Btz] [ Bzz] −1 ( [ Czt] −j k0 n0 [ Bzt])+ ( [ Ctz] +j k0 n0 [ Btz]) [ Bzz] −1 [ Bzt]
(17)
[ A12]= [ Btt] − [ Btz] [ Bzz] −1 [ Bzt]
(18)
[ B11]= − [ Att] − k02 n02 [ Btt] − (j [ Ctz] − k0 n0 [ Btz]) [ Bzz] −1 (j [ Czt] + k0 n0 [ Bzt]).
(19)

The sparsity of the matrices from Eqs. (5)–(11) can be employed in the calculation of these matrices. Discretization of the derivative in Eq. (14) yields following recurrence scheme for propagation of the transverse field in the z-direction:

[Y]i [ ϕt d ϕtdz] i+1= [Z]i [ ϕt d ϕtdz]i
(20)

with i the iteration number and

[Y]i= [A]i+ϑΔz [B]i
(21)
[ Z]i= [A]i − ( 1−ϑ)Δz [B]i,
(22)

where Δz is the propagation step and ϑ is a parameter that controls the stability of the scheme. The subscript i for the system matrices [Y] and [Z] in Eq. (20) can be omitted for structures with a cross section which is invariant in the z-direction (e.g. waveguides).

In contrast to other reported finite element BPMalgorithms that assume a diagonal permittivity tensor [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

,6

D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, “Mixed finite element beam propagation method,” J. Lightwave Technol. 16, 1336–1342 (1998). [CrossRef]

], the recurrence scheme presented in Eq. (20) is applicable to general anisotropic materials because all elements of the dielectric tensor are included in the calculation of [A] i and [B] i . The longitudinal field component has been eliminated in Eq. (12) to obtain the recurrence scheme for the transverse field in Eq. (20). Koshiba et al. [11

M. Koshiba and K. Inoue, “Simple and efficient finite-element analysis of microwave and optical waveguides,” IEEE Trans. Microwave Theory Tech. 40, 371–377 (1992). [CrossRef]

] have used a similar elimination to obtain an eigenvalue problem for the calculation of waveguide eigenmodes in which a diagonal permittivity tensor is assumed. The loss of the sparsity of the matrices in Eqs. (17)–(19) is the price paid to obtain such a recurrence scheme which includes the full dielectric tensor and avoids the use of a first order approximation of the wave equation to make the method more accurate.

Including not only the transverse field but its derivative with respect to z too in the propagation vector avoids the use of a Padé [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

, 6

D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, “Mixed finite element beam propagation method,” J. Lightwave Technol. 16, 1336–1342 (1998). [CrossRef]

, 12

G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. 17, 1743–1745 (1992). [CrossRef] [PubMed]

] approximation or some other first order approximation to deal with the second order derivatives that are present in Eq. (13). Instead, the initial derivative of the input transverse field with respect to z is required at the point z=0 to start recursive application of Eq. (20) to calculate the evolution of the electric field upon propagation. In order to benefit from the higher accuracy of the second order system compared to first order approximations of the wave equation, it is important to calculate this initial derivative of the transverse field accurately.

First, the input transverse field is discretized on a mesh perpendicular to the propagation direction. A preliminary paraxial propagation of this transverse field is considered over a very short distance Δz′. In the paraxial approximation, the second order derivative in Eq. (13) is neglected which leads to a differential equation similar to Eq. (14) (with ϕt as the only variable). The fact that this preliminary propagation is paraxial does not restrict the validity of further wide angle propagation as this step is only performed to calculate a start value of the derivative of the transverse field at z=0, which will be further refined before the actual wide angle propagation is started. A paraxial recurrence scheme similar to Eq. (20) can be obtained by replacing [Y] i and [Z] i by their paraxial equivalents [Y′] and [Z′] respectively:

[ Y′]= [ A11]+ϑΔz [ B11]
(23)
[ Z′]= [ A11] − ( 1−ϑ)Δz [ B11].
(24)

Using the transverse field at z=z′ obtained from the paraxial propagation step, a start value for the derivative of the input field at z=0 can be calculated as a forward difference:

d ϕtdz∣ z=0 ≈ ϕt ∣ z=Δ z′ − ϕt∣ z=0 Δ z′,
(25)

if Δz′ is chosen to be sufficiently small. Next, this start value for the derivative at z=0 is used together with the input transverse field in the wide angle recurrence scheme from Eq. (20) to consider the propagation of this beam over a distance Δz. The transverse field obtained at z=Δz is then again used to calculate a next, more accurate, value for the derivative at z=0 similar as in Eq. (25). It can be shown that if this wide angle iteration procedure is repeated for a sufficient number N of steps, the initial derivative at z=0 and the transverse field at z=Δz will converge to consistent values. The actual wide angle beam propagation can then be started with these values, taking advantage of the higher accuracy of the second order system.

The evolution of the electric field upon (wide angle) propagation of the input optical field through the structure of interest is calculated from iterative application of Eq. (20). This yields the transverse field ϕt and its derivative with respect to z in subsequent planes separated by Δz. The longitudinal field component ϕz can be calculated from Eq. (12). To realize high stability of the propagation algorithm over large propagation distances, it is possible (but not required) to carry out a preliminary wide angle propagation step between positions z 0 and z 0+Δz to calculate the central derivative of the transverse field with respect to z at z 0. This derivative can then be inserted into the propagation vector at z 0 before the wide angle propagation step is carried out towards z 0+Δz.

2.3. Implementation

The presented finite element BPM was implemented in Matlab and takes advantage of the efficient operations on sparse matrices. The mesh file is generated by a dedicated program such as GiD [13

GiD, the personal pre and post processor, http://gid.cimne.upc.es/.

] and the mesh is imported into the Matlab workspace. The submatrices in Eqs. (5)–(11) are calculated using the Gauss Quadrature Rule for numerical integration [9

J. Jin, The finite element method in electromagnetics , 2nd edition (Wiley, New York US, 2002).

].

Running the BPM on a personal computer with 4GB of Random Access Memory allows to simulate structures consisting of about 500 LT/QN elements of which the edge dimensions are about the order of the light wavelength. Typical parameter values used in our BPM simulations are: the propagation step Δz≈λ, the PML thickness di =1.5µm (i=1 to 4), the stability parameter ϑ=1 and maximum attenuation αmax =12.5 in the PML regions [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

]. The total calculation time to simulate light propagation in general or inhomogeneous anisotropic materials (full dielectric tensor) over e.g. 100µm for a mesh consisting of 500 LT/QN elements is typically 10 minutes on a 2.5GHz Intel Core 2 Duo CPU.

A compiled evaluation version of the BPM program will be made available for download at www.elis.ugent.be/ELISgroups/lcd/research/bpm.php.

3. Applications

In this Section practical applications of the wide angle BPM algorithm are presented to demonstrate the wide applicability of the method. First, the propagation of a laser beam in a homogeneous uniaxial material is discussed and the simulation results are compared with theoretical expectations and the numerical error is evaluated. Next, the tunability of the coupling characteristics of a directional coupler based on liquid crystals is evaluated. The liquid crystal orientation is obtained from simulation with a finite element solver [14

R. James, E. Willman, F. A. Fernandez, and S. E. Day, “Finite-Element Modeling of Liquid-Crystal Hydrodynamics With a Variable Degree of Order,” IEEE Trans. Electron Devices 53, 1575–1582 (2006). [CrossRef]

] and the resulting director profile is considered in the BPM analysis of the coupler. This example illustrates how the BPM can be used for optical analysis of liquid crystal devices, taking into account the elastic distortion of the liquid crystal e.g. due to fringing field effects. A third example shows the diffraction of a plane wave propagating through structures containing micrometer variations in permittivity, such as the pixels in a spatial light modulator (SLM). Unless otherwise stated, the simulation parameters are set as described in Subsection 2.3.

3.1. Propagation of a laser beam in a uniaxial layer

In this example the propagation of a gaussian laser beam in a homogeneous uniaxial layer is simulated using the BPM algorithm. Considering a homogeneous layer allows comparison of simulated results and theoretical expectations in order to verify the accuracy of the proposed model. Firstly, it is shown that polarization rotation is induced by propagation through a transversely anisotropic material. Next, the difference between the direction of the flow of optical energy (Poynting vector) and the wave vector is shown upon propagation in a material with longitudinal dielectric anisotropy.

3.1.1. Polarization rotation

The polarization rotation of an electromagnetic plane wave upon propagation through an anisotropic optical material is well-known. As a first illustration of the use of the developed BPM algorithm, the polarization rotation of a linearly polarized beam with a gaussian transverse intensity profile at z=0 is calculated upon propagation through a homogeneous uniaxial layer. The Ex and Ey field components of the input optical field with wavelength λ=1µm are in phase at z=0, leading to a linear polarization state parallel to the first bisector as illustrated in Figs. 2(a)–2(c). The Ez component of the input optical field shown in Fig. 2(d) is calculated from Eq. (12) where the initial derivative of the transverse field dϕtdzat z=0 at z=0 is obtained after a paraxial propagation of the transverse field over a short distance Δz′=0.1mm and N=10 wide angle iterations, as explained in Subsection 2.2.

The input beam shown in Fig. 2 at z=0 propagates through a homogeneous uniaxialmaterial with refractive indices nx =1.5, ny =1.6 and nz =1.5. There is a weak broadening of the gaussian intensity profile due to diffraction, and as a result the amplitude of the transverse field components Ex and Ey decrease slowly during propagation. Furthermore, Ex and Ey will propagate at different phase velocities c/nx resp. c/ny through the medium, which leads to a phase delay Γ=2π/λ(ny−nx )d between these components after propagation over a distance d. The change in the phase delay G will alter the polarization state of the light. This is illustrated in Figs. 3(a)–3(d) which show the evolution of the Ez component upon propagation through the uniaxial layer. The BPM is run with a propagation step of Δz=0.25µm on a mesh comprising 384 elements.

Fig. 2. Input optical field at z=0µm: (a) transverse electric field, (b) Ex component, (c) Ey component, (d) Ez component.
Fig. 3. Evolution of Ez upon propagation through the uniaxial layer: (a) d=1.25µm, (b) d=2.5µm, (c) d=3.75µm, (d) d=5mm.

For linearly polarized light there is a symmetry axis through the lobes of the Ez field component parallel to the polarization state as can be seen by comparing Figs. 2(a) and 2(d). After propagation over a distance d=1.25µm, the phase difference between Ex and Ey becomes Γ=π/4 which leads to an elliptical polarization state with long axis parallel to the first bisector of the transverse plane, which can be seen from the intensity profile of Ez in Fig. 3(a). For d=2.5µm, the Ex and Ey components are in quadrature with Γ=π/2, which leads to a circular polarization (see Fig. 3(b)). After d=3.75µm the phase difference Γ=3π/4 leads to an elliptic polarization state with long axis parallel to the −45° bisector (Fig. 3(c)). Finally, after propagation over a distance d=5µm, Ex and Ey become opposite in sign and a linear polarization state is obtained which is perpendicular to the initial polarization (Fig. 3(d)). It is clear that the simulated evolution of the field profile obtained with the BPM agrees well with theoretical expectations.

3.1.2. Deviation angle d between the Poynting vector and the wave vector

It can be shown [15

J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect,” Opt. Quantum Electron. 37, 95–106 (2005). [CrossRef]

] that the direction of the flow of optical energy (Poynting vector) deviates from the wave vector in the material if dielectric anisotropy is present in the longitudinal direction (e.g. εzy ≠0). For example, consider a plane wave with wave vector k̄i at the interface between air and an anisotropic material as shown in Fig. 4(a). In this example, the superscripts i and t denote quantities related to the incident and transmitted plane wave, respectively.

If k̄i is oriented along the z-axis, the wave vector of the transmitted plane wave k̄t will only have a z-component. According to the boundary conditions from Maxwell’s equations, continuity of the normal component of the displacement field D̄ t =ε̿Ēt at the interface is required. If the electric field of the incident plane wave is oriented along the y-axis (Ei y ≠0,Ei x =Ei z =0), this yields:

0= Dzt
(26)
= εzy Eyt+ εzz Ezt.
(27)

The magnetic field H̄t of the plane wave transmitted to the anisotropic medium is oriented parallel to the x-axis (Ht x ≠0,Ht y =Ht z =0). Equation 27 shows that the electric field Ēt of the transmitted wave will no longer be perpendicular to k̄t . As a consequence, the Poynting vector S̄ t =Ēt ×H̄t , representing the direction of the flow of optical energy, will not be along the z-axis. The angle d between the Poynting vector S̄ and the wave vector k̄t shown in Fig. 4(a) is calculated as:

tanδ= Syt Szt = − Ezt Hxt Eyt Hxt = εzy εzz.
(28)

To illustrate the occurrence of the deviation angle δ, we calculated the propagation of an optical field with wave vector parallel to the z-axis in a homogeneous uniaxial material with ordinary and extra-ordinary refractive indices n o=1.5 and ne =1.6, respectively. The input field has a gaussian intensity profile for Ey , linear vertical polarization (Ex =0) and wavelength λ=1µm. The inclination angle of the crystalline c-axis of the uniaxial layer relative to the y-axis (as indicated in Fig. 4(b)) is θ=π/4. The azimuth of the c-axis relative to the x-axis in the xz-plane (Fig. 4(b)) is ϕ=π/2. The dielectric tensor (at optical frequencies) for given inclination θ and azimuth ϕ is:

ε̿= [ ε⊥+Δε sin2θ cos2ϕ Δεsinθcosθcosϕ Δε sin2θcosϕsinϕ Δεsinθcosθcosϕ ε⊥+Δε cos2θ Δεsinθcosθsinϕ Δε sin2θcosϕsinϕ Δεsinθcosθsinϕ ε⊥+Δε sin2θ sin2ϕ],
(29)

with ε ⊥=n 2 0, ε ‖=n 2 e and Δε=ε ‖−ε ⊥. The dielectric components εzy and εzz can be evaluated with Eq. (29) and the corresponding deviation angle δ according to Eq. (28) is given by:

δ=arctan ( εzy εzz) ≈3.688°.
(30)

Figure 5 shows the intensity profile of the optical field (a) before and (b) after propagation over 10µm through the uniaxial layer, simulated with the BPM for a propagation step Δz=0.25µm.

Fig. 4. (a) Wave vectors of incident and transmitted plane waves at an interface between air and an anisotropic material ε̿. (b) Definition of inclination θ and azimuth ϕ angles.
Fig. 5. Intensity profile of the gaussian optical field (a) at z=0µm and (b) at z=10µm.

Alongside a broadening of the optical field due to diffraction, comparing Figs. 5(a) and 5(b) reveals a vertical displacement of the beam upon propagation. The simulated displacement Δy=0.65µm obtained for a mesh consisting of 248 elements is in good agreement with the predicted deviation angle (calculated for plane waves) δ=3.688° which implies a vertical translation over 0.645µm after propagation over a distance of 10µm. The accuracy can be further improved by increasing the number of elements in the mesh. Figure 6 shows the relative error on the beam displacement from the theoretical value after propagation over a distance of 10µm as a function of the number of elements used for a square 10µm by 10µm mesh.

Fig. 6. The relative error on the beam displacement after propagation over 10µm as a function of the number of elements.

As expected, the numerical error on the beam displacement decreases with the number of elements in Fig. 6. This error is about 10−5 when 300 elements are used and below 10−6 for 500 elements. The numerical accuracy of the perfectly matched layers that are applied is evaluated in [5

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

]. The value of the maximum attenuation in the PML regions as defined in Subsection 2.3 is set according to this analysis to have a small numerical error due to truncation of the computational window. The good agreement between the simulated results and theoretical expectations shows that the light propagation in anisotropic media is modeled accurately using the presented BPM algorithm.

3.2. A tunable directional coupler based on liquid crystals

If two waveguides are sufficiently close to each other, light can be coupled fromone to the other. The interchange of optical power between the two waveguides can be used to make optical couplers and switches [16

J. C. Campbell, F. A. Blum, D. W. Shaw, and K. L. Lawlay, “GaAs Electro-optic directional coupler switch,” Appl. Phys. Lett. 27, 202–205 (1975). [CrossRef]

]. In this Subsection we consider a directional coupler consisting of two square waveguides and a transverse structure cross section as in Fig. 1(a) with dimensions l 1=2µm, l 2=3µm, l 3=4µm, l 4=11µm. The two waveguides each have an area of 1mm2 and they are separated by a distance d=1µm. To realize a tunable coupler, the waveguides with refractive index n 1=1.7 are surrounded by a layer of uniaxial liquid crystal with ordinary and extraordinary refractive index no =1.5 and ne =1.6, respectively. The substrates have a refractive index n 3=1.5. The orientation of the liquid crystal molecules can be controlled because of their dielectric anisotropy by applying a voltage between the electrodes E1 and E2 indicated in Fig. 1(a). A homogeneous orientation of the liquid crystal layer in the 0V state can be realized by using an appropriate alignment layer at the interfaces of the liquid crystal layer with the surrounding media ε 1 and ε 3. In this example, the liquid crystal molecules have an initial planar alignment with inclination (pretilt) θ=88° and azimuth (pretwist) ϕ=90°, where θ and ϕ are defined as in Fig. 4(b).

The orientation of the liquid crystal is found by minimizing the Landau-de-Gennes (LdG) free energy functional [17

P. G. de Gennes and J. Prost, The Physics of Liquid Crystals , 2nd edition (Clarendon, Oxford UK, 1993).

] that comprises bulk, elastic, surface and electrostatic contributions. This approach allows for variations in the order parameter of the liquid crystal, required for the proper description of disclinations. To calculate the orientation of the liquid crystal in the coupler, a two-dimensional finite element implementation as described in [14

R. James, E. Willman, F. A. Fernandez, and S. E. Day, “Finite-Element Modeling of Liquid-Crystal Hydrodynamics With a Variable Degree of Order,” IEEE Trans. Electron Devices 53, 1575–1582 (2006). [CrossRef]

] has been used. Figures 7 and 8 show the calculated director profiles in the directional coupler when voltages of 0V and 7V are applied between the electrodes.

As can be seen from Figs. 7 and 8, the initial planar alignment of the liquid crystal can be changed to a nearly vertical alignment by applying a potential difference of 7V across the electrodes. Because of the optical anisotropy of the liquid crystal, changing the orientation of the molecules alters the optical properties of the liquid crystal surrounding the waveguides of the directional coupler. In each element of the meshed liquid crystal layer, the components of the dielectric tensor (at optical frequencies) can be calculated from the tilt θe and twist ϕe angles within the element by using Eq. (29).

As the coupling characteristics of the directional coupler are determined by the effective indices of the supermodes propagating in the waveguides, changing the optical properties of the surrounding liquid crystal layer will alter the coupling behavior. Therefore, the presented structure will behave as a voltage controllable tunable directional coupler.

Fig. 7. Director profile of the liquid crystal molecules when 0V is applied. The tilt angle θ of the director in the liquid crystal layer is indicated in the color bar.
Fig. 8. Director profile of the liquid crystal molecules when 7V is applied. The tilt angle θ of the director in the liquid crystal layer is indicated in the color bar.

To simulate the light propagation through the coupler, we apply the presented finite element BPM. It follows from Eq. (29) that it is indeed required to incorporate the full dielectric tensor in each point independently to correctly simulate the behavior of the inhomogeneous liquid crystal layer because in each element, depending on the values of θe and ϕe obtained from the director simulation, all nine elements in this tensor can be non-zero. Therefore, an arbitrary orientation of the liquid crystal directors is allowed here, in contrast to the method proposed in [4

Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. 23, 2014–2019 (2006). [CrossRef]

] which requires all directors to be aligned in the xy-plane. Compared to the finite difference methods reported in [3

A. J. Davidson and S. J. Elston, “Three-dimensional beam propagation model for the optical path of light through a nematic liquid crystal,” J. Mod. Opt. 53, 979–989 (2006). [CrossRef]

,4

Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. 23, 2014–2019 (2006). [CrossRef]

,18

N. Amarasinghe, E. Gartland, and J. Kelly, “Modeling optical properties of liquid-crystal devices by numerical solution of time-harmonic Maxwell equations,” J. Opt. Soc. Am. 21, 1344–1361 (2004). [CrossRef]

] for the study of light propagation in liquid crystal devices, the use of a finite element implementation allows using irregular meshes of arbitrary shape which can be optimized to model the structure in an efficient way by concentrating elements in the areas of interest.

To determine the coupling characteristics of the directional coupler, the propagation of an input optical field [0,Ey,Ez ] with a vertical linear polarization and wavelength λ=1.5µm was simulated in the z-direction for a propagation step Δz=1µm. The Ey field component has a gaussian intensity profile centered in the left waveguide and does not overlap with the second waveguide. Figures 9 and 10 show the evolution of the Ey and Ez components upon propagation through the directional coupler when 0V and 7V, respectively, is applied over the electrodes. The computational window shown in Fig. 1(a) is divided into 770 triangular finite elements.

Figures 9 and 10 show the coupling of the input optical field from the left waveguide to the right waveguide. When no voltage is applied (Fig. 9), the optical power in both waveguides is about equal after a propagation distance d=50µm. After propagation over 105µm all optical power is transferred to the right waveguide. When 7V is applied across the electrodes, waveguide coupling occurs similarly to the 0V case. However, comparison of Figs. 9 and 10 shows that the coupling length is decreased because of the reduced inclination of the liquid crystalline molecules (see Fig. 8). Indeed, after a propagation distance of 50µm already more than 50% of the optical power is coupled to the right waveguide. After propagation over 90µm all optical power has transferred from one waveguide to the other. This illustrates the tunability of the presented directional coupler: by applying 7V over the liquid crystal layer, the coupling length can be decreased from 105µm to 90µm. Further optimization of the geometry of the structure and material parameters may increase the tunability of the structure. This example illustrates how the presented BPM algorithm can be combined with a liquid crystal director simulation.

Fig. 9. Evolution of the electric field components Ey (top) and Ez (bottom) upon propagation through the directional coupler when no voltage is applied. The contour plots are normalized to have equal maximum values for each frame.
Fig. 10. Evolution of the electric field components Ey (top) and Ez (bottom) upon propagation through the directional coupler when 7V is applied over the electrodes. The contour plots are normalized to have equal maximum values for each frame.

3.3. Diffraction in a spatial light modulator (SLM)

As illustrated in the previous example, the developed BPM is an interesting tool to simulate the light propagation in liquid crystal devices. As well as being able to model continuous layers of liquid crystalline materials, such as the directional coupler possesses, the BPM can also be used to simulate light propagation through structures consisting of pixels as in liquid crystal displays (LCDs) or spatial light modulators (SLMs). Applying the described BPM becomes relevant for the simulation of structures in which the pixel dimensions approach the wavelength of light. In this case a rigorous vector method is required to accurately model light diffraction, which plays a crucial role. Recently, small pixel SLMs are gaining interest for holographic display applications [19

E. Buckley, “Holographic Laser Projection Technology,” SID Int. Symp. Digest Tech. Papers 39, 1074–1079 (2008). [CrossRef]

] as phase modulating components. Application of the presented BPM can provide detailed information on the near-field profile of light after propagation.

As an example, we simulate the propagation in the z-direction of a plane wave [0,Ey ,0] with vertical linear polarization and wavelength λ=1µm through two structures each consisting of 4 pixels as shown in Fig. 11 with dimensions lp =4µm and lp =2µm, respectively. To illustrate the influence of the pixel size only on light diffraction whilst excluding boundary effects due to the elastic distortion of the liquid crystal between neighboring pixels, ideal pixels are considered which are modeled as homogeneous uniaxial materials with refractive indices no =1.5 and ne =1.6. For the top left and bottom right pixel, which are in the ON state, the crystalline c-axis is oriented along the first bisector of the xy-plane (corresponding to θ=π/4 and ϕ=0 in terms of Eq. (29)). In the top right and bottom left pixel, which are in the OFF state, the c-axis is parallel to the z-axis (θ=π/2 and ϕ=π/2).

Fig. 11. Sketch of the four pixel structure, indicating the ON and OFF pixels.

If we consider the propagation of the plane wave [0,Ey ,0] through the structure shown in Fig. 11, a polarization change will occur in the ON pixels in an analogous way as described in Subsection 3.1.1 while the polarization state will remain unchanged in the OFF pixels. Fig. 12 shows the amplitude of the near-field Ex component obtained by simulating the propagation of the Ey plane wave through the structure over d=5µm for pixel dimensions (a) 4µm by 4µm and (b) 2µm by 2µm. The simulation results were obtained for a propagation step Δz=0.25µm and the simulation window was divided into 432 triangular elements.

The contour plots of Ex in Fig. 12 show that the original vertical linear polarization of the input plane wave (Ex =0) is indeed rotated over 90° after propagation through the ON pixels over 5µm to obtain a horizontal linear polarization. As expected, the initial vertical polarization largely remains unchanged in the OFF pixels (Ex ≈0). Comparing Figs. 12(a) and 12(b) illustrates that the pixel dimensions have an important influence on the near-field diffraction of the light after propagation through the structure. For the 4µm by 4µm pixels, the resulting amplitude profile of Ex consists of several peaks due to elementary Fresnel diffraction of a plane wave. Due to the smaller pixel dimensions, the Ex component has a more confined diffraction profile in the 2µm by 2µm pixels. It can be seen from Fig. 12 that due to light diffraction of Ex in the ON pixels, this component spreads out to partially cover the OFF pixels. This effect is already present for the 4µm by 4µm pixels (Fig. 12(a)) and degrades the contrast of the 2µm by 2µm pixels to a large extent (Fig. 12(b)). This fundamental degradation in contrast puts a limit on the pixel size for high quality amplitude modulation SLMs.

Fig. 12. Amplitude of the near-field Ex component obtained after propagation of a plane wave [0,Ey ,0] over 5µm through the 4-pixel structure for (a) 4µm by 4µm and (b) 2µm by 2µm pixels.

It is clear from this example that the presented BPM can provide detailed information on the near-field pattern to aid the design of advanced structures consisting of several pixels such as LCDs or SLMs.

4. Conclusions

A full-vectorial, finite element wide angle beam propagation method has been presented that uses higher order vector elements and is applicable to general anisotropic dielectric materials. As no approximations are made concerning the shape and elements of the dielectric tensor, the method is well-suited to simulate light propagation in inhomogeneous anisotropic materials such as liquid crystals. Employing a finite element scheme offers the flexibility to model arbitrarily curved structures. The general applicability of the presented method is illustrated by four examples showing the propagation of a laser beam in homogeneous uniaxial materials with transverse and longitudinal dielectric anisotropy, the tunability of a liquid crystal based directional coupler and the near-field diffraction of a plane wave in a structure containing micrometer-scale variations similar to the pixels of a spatial light modulator. The combination of the efficiency of the beam propagation algorithm, the flexibility of the finite element implementation and the wide applicability to general anisotropic materials makes the presented method an interesting tool for the optical analysis of photonic components.

Acknowledgments

Pieter Vanbrabant is a PhD Fellow of the Research Foundation - Flanders (FWO Vlaanderen) and Jeroen Beeckman is a Postdoctoral Fellow of the same institution. The work has been carried out in the framework of the Interuniversity Attraction Pole (IAP) project Photonics@be of the Belgian Science Foundation.

References and links

1.

J. M. López-Doña, J. G. Wangüemert-Pérez, and I. Molina-Fernández, “Fast-fourier-based three-dimensional full-vectorial beam propagation method,” IEEE Photonics Technol. Lett. 17, 2319–2321 (2005). [CrossRef]

2.

C. Ma and E. Van Keuren, “A three-dimensional wide-angle BPM for optical waveguide structures,” Opt. Express 15, 402–407 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?uri=OE-15-2-402. [CrossRef] [PubMed]

3.

A. J. Davidson and S. J. Elston, “Three-dimensional beam propagation model for the optical path of light through a nematic liquid crystal,” J. Mod. Opt. 53, 979–989 (2006). [CrossRef]

4.

Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. 23, 2014–2019 (2006). [CrossRef]

5.

K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]

6.

D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, “Mixed finite element beam propagation method,” J. Lightwave Technol. 16, 1336–1342 (1998). [CrossRef]

7.

M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). [CrossRef]

8.

F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microw. Guid. Wave Lett. 8, 223–225 (1998). [CrossRef]

9.

J. Jin, The finite element method in electromagnetics , 2nd edition (Wiley, New York US, 2002).

10.

J. Beeckman, R. James, F. A. Fernandez, W. De Cort, P. J. M. Vanbrabant, and K. Neyts, “Calculation of fully anisotropic liquid crystal waveguide modes,” accepted for publication in J. Lightwave Technol. (2009). [CrossRef]

11.

M. Koshiba and K. Inoue, “Simple and efficient finite-element analysis of microwave and optical waveguides,” IEEE Trans. Microwave Theory Tech. 40, 371–377 (1992). [CrossRef]

12.

G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. 17, 1743–1745 (1992). [CrossRef] [PubMed]

13.

GiD, the personal pre and post processor, http://gid.cimne.upc.es/.

14.

R. James, E. Willman, F. A. Fernandez, and S. E. Day, “Finite-Element Modeling of Liquid-Crystal Hydrodynamics With a Variable Degree of Order,” IEEE Trans. Electron Devices 53, 1575–1582 (2006). [CrossRef]

15.

J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect,” Opt. Quantum Electron. 37, 95–106 (2005). [CrossRef]

16.

J. C. Campbell, F. A. Blum, D. W. Shaw, and K. L. Lawlay, “GaAs Electro-optic directional coupler switch,” Appl. Phys. Lett. 27, 202–205 (1975). [CrossRef]

17.

P. G. de Gennes and J. Prost, The Physics of Liquid Crystals , 2nd edition (Clarendon, Oxford UK, 1993).

18.

N. Amarasinghe, E. Gartland, and J. Kelly, “Modeling optical properties of liquid-crystal devices by numerical solution of time-harmonic Maxwell equations,” J. Opt. Soc. Am. 21, 1344–1361 (2004). [CrossRef]

19.

E. Buckley, “Holographic Laser Projection Technology,” SID Int. Symp. Digest Tech. Papers 39, 1074–1079 (2008). [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(130.2790) Integrated optics : Guided waves
(350.5500) Other areas of optics : Propagation

ToC Category:
Optical Devices

History
Original Manuscript: February 10, 2009
Revised Manuscript: April 30, 2009
Manuscript Accepted: May 8, 2009
Published: June 16, 2009

Citation
Pieter J. M. Vanbrabant, Jeroen Beeckman, Kristiaan Neyts, Richard James, and F. Anibal Fernandez, "A finite element beam propagation method for simulation of liquid crystal devices," Opt. Express 17, 10895-10909 (2009)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-13-10895


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. J. M. Lopez-Dona, J. G. Wanguemert-Perez, and I. Molina-Fernandez, "Fast-fourier-based three-dimensional full-vectorial beam propagation method," IEEE Photonics Technol. Lett. 17, 2319-2321 (2005). [CrossRef]
  2. C. Ma and E. Van Keuren, "A three-dimensional wide-angle BPM for optical waveguide structures," Opt. Express 15, 402-407 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?uri= OE-15-2-402. [CrossRef] [PubMed]
  3. A. J. Davidson and S. J. Elston, "Three-dimensional beam propagation model for the optical path of light through a nematic liquid crystal," J. Mod. Opt. 53, 979-989 (2006). [CrossRef]
  4. Q. Wang, G. Farrell, and Y. Semenova, "Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method," J. Opt. Soc. Am. 23, 2014-2019 (2006). [CrossRef]
  5. K. Saitoh and M. Koshiba, "Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides," J. Lightwave Technol. 19, 405-413 (2001). [CrossRef]
  6. D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, "Mixed finite element beam propagation method," J. Lightwave Technol. 16, 1336-1342 (1998). [CrossRef]
  7. M. Koshiba and Y. Tsuji, "Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems," J. Lightwave Technol. 18, 737-743 (2000). [CrossRef]
  8. F. L. Teixeira and W. C. Chew, "General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media," IEEE Microw. Guid. Wave Lett. 8, 223-225 (1998). [CrossRef]
  9. J. Jin, The finite element method in electromagnetics, 2nd edition (Wiley, New York US, 2002).
  10. J. Beeckman, R. James, F. A. Fernandez, W. De Cort, P. J. M. Vanbrabant, and K. Neyts, "Calculation of fully anisotropic liquid crystal waveguide modes," accepted for publication in J. Lightwave Technol. (2009). [CrossRef]
  11. M. Koshiba and K. Inoue, "Simple and efficient finite-element analysis of microwave and optical waveguides," IEEE Trans. Microwave Theory Tech. 40, 371-377 (1992). [CrossRef]
  12. G. R. Hadley, "Multistep method for wide-angle beam propagation," Opt. Lett. 17, 1743-1745 (1992). [CrossRef] [PubMed]
  13. GiD, the personal pre and post processor, http://gid.cimne.upc.es/.
  14. R. James, E. Willman, F. A. Fernandez, and S. E. Day, "Finite-Element Modeling of Liquid-Crystal Hydrodynamics With a Variable Degree of Order," IEEE Trans. Electron Devices 53, 1575-1582 (2006). [CrossRef]
  15. J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, "Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect," Opt. Quantum Electron. 37, 95-106 (2005). [CrossRef]
  16. J. C. Campbell, F. A. Blum, D. W. Shaw, and K. L. Lawlay, "GaAs Electro-optic directional coupler switch," Appl. Phys. Lett. 27, 202-205 (1975). [CrossRef]
  17. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, 2nd edition (Clarendon, Oxford UK, 1993).
  18. N. Amarasinghe, E. Gartland, and J. Kelly, "Modeling optical properties of liquid-crystal devices by numerical solution of time-harmonic Maxwell equations," J. Opt. Soc. Am. 21, 1344-1361 (2004). [CrossRef]
  19. E. Buckley, "Holographic Laser Projection Technology," SID Int. Symp. Digest Tech. Papers 39, 1074-1079 (2008). [CrossRef]

Cited By

 Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited