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

Optics Express, Vol. 17, Issue 13, pp. 10895-10909 (2009)

http://dx.doi.org/10.1364/OE.17.010895

Acrobat PDF (2928 KB)

### 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

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]

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

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

### 2.1. Finite element discretization of the wave equation

*E*̄(

*x,y, z*), as derived from the Maxwell equations:

*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

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]

*µ*̿

^{−1}and

*ε*̿ with

*p*̿ and

*q*̿ respectively. These are permittivities and permeabilities modified to include the absorption coefficients in the eight PML regions [5

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

*E*̄ 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̄

*and a phase factor exp(−*

_{z}*jk*

_{0}

*n*

_{0}

*z*), where

*z*is the direction of light propagation:

*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:

*ϕ*and

^{e}_{t}*ϕ*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)

^{e}_{z}*N*and

_{x}*N*and quadratic normal (QN)

_{y}*L*shape functions for interpolation of the transverse and longitudinal field [7

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]

*ε*̿ is maintained and a Galerkin procedure [9] is applied to the wave equation from Eq. (1). Substituting Eqs. (2) and (3) into Eq. (1) leads to the beam propagation expression:

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]

*∂p*is sufficiently small, the submatrices in Eq. (4) can be evaluated as:

_{ii}/∂z*ϕ*into the first row of Eq. (4) yields a second order differential equation for the transverse field component:

_{z}### 2.2. Wide angle beam propagation algorithm

*z*-direction, Eq. (13) is formally rewritten in the form of a first order differential equation:

*A*and

*B*are defined as:

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]

*A*

_{11}], [

*A*

_{12}] and [

*B*

_{11}] which are defined as:

*z*-direction:

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

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]

*A*]

*and [*

_{i}*B*]

*. 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*

_{i}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]

*z*too in the propagation vector avoids the use of a Padé [5

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]

*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.

*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

*ϕ*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

_{t}*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*]

*and [*

_{i}*Z*]

*by their paraxial equivalents [*

_{i}*Y*′] and [

*Z*′] respectively:

*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:

*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.

*ϕ*and its derivative with respect to

_{t}*z*in subsequent planes separated by Δ

*z*. The longitudinal field component

*ϕ*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}*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

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

*z*≈

*λ*, the PML thickness

*d*=1.5

_{i}*µ*m (

*i*=1 to 4), the stability parameter ϑ=1 and maximum attenuation

*α*=12.5 in the PML regions [5

_{max}5. K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. **19**, 405–413 (2001).
[CrossRef]

*µ*m for a mesh consisting of 500 LT/QN elements is typically 10 minutes on a 2.5GHz Intel Core 2 Duo CPU.

## 3. Applications

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]

### 3.1. Propagation of a laser beam in a uniaxial layer

#### 3.1.1. Polarization rotation

*z*=0 is calculated upon propagation through a homogeneous uniaxial layer. The

*E*and

_{x}*E*field components of the input optical field with wavelength

_{y}*λ*=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

*E*component of the input optical field shown in Fig. 2(d) is calculated from Eq. (12) where the initial derivative of the transverse field

_{z}*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.

*z*=0 propagates through a homogeneous uniaxialmaterial with refractive indices

*n*=1.5,

_{x}*n*=1.6 and

_{y}*n*=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

_{z}*E*and

_{x}*E*decrease slowly during propagation. Furthermore,

_{y}*E*and

_{x}*E*will propagate at different phase velocities

_{y}*c/n*resp.

_{x}*c/n*through the medium, which leads to a phase delay Γ=2

_{y}*π/λ*(

*n*)

_{y}−n_{x}*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

*E*component upon propagation through the uniaxial layer. The BPM is run with a propagation step of Δ

_{z}*z*=0.25

*µ*m on a mesh comprising 384 elements.

*E*field component parallel to the polarization state as can be seen by comparing Figs. 2(a) and 2(d). After propagation over a distance

_{z}*d*=1.25

*µ*m, the phase difference between

*E*and

_{x}*E*becomes

_{y}*Γ*=π/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

*E*in Fig. 3(a). For

_{z}*d*=2.5

*µ*m, the

*E*and

_{x}*E*components are in quadrature with

_{y}*Γ*=π/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,

*E*and

_{x}*E*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.

_{y}#### 3.1.2. Deviation angle d between the Poynting vector and the wave vector

*k̄*is oriented along the

^{i}*z*-axis, the wave vector of the transmitted plane wave

*k̄*will only have a

^{t}*z*-component. According to the boundary conditions from Maxwell’s equations, continuity of the normal component of the displacement field

*D*̄

*=*

^{t}*ε̿Ē*at the interface is required. If the electric field of the incident plane wave is oriented along the

^{t}*y*-axis (

*E*≠0,

^{i}_{y}*E*=

^{i}_{x}*E*=0), this yields:

^{i}_{z}*H̄*of the plane wave transmitted to the anisotropic medium is oriented parallel to the

^{t}*x*-axis (

*H*≠0,

^{t}_{x}*H*=

^{t}_{y}*H*=0). Equation 27 shows that the electric field

^{t}_{z}*Ē*of the transmitted wave will no longer be perpendicular to

^{t}*k̄*. As a consequence, the Poynting vector

_{t}*S*̄

*=*

^{t}*Ē*×

^{t}*H̄*, representing the direction of the flow of optical energy, will not be along the

^{t}*z*-axis. The angle d between the Poynting vector

*S*̄ and the wave vector

*k̄*shown in Fig. 4(a) is calculated as:

_{t}*δ*, 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

*n*=1.6, respectively. The input field has a gaussian intensity profile for

_{e}*E*, linear vertical polarization (

_{y}*E*=0) and wavelength

_{x}*λ*=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:

*ε*

_{⊥}=

*n*

^{2}

_{0},

*ε*

_{‖}=

*n*

^{2}

*and Δ*

_{e}*ε*=

*ε*

_{‖}−

*ε*

_{⊥}. The dielectric components

*ε*and

_{zy}*ε*can be evaluated with Eq. (29) and the corresponding deviation angle

_{zz}*δ*according to Eq. (28) is given by:

*z*=0.25

*µ*m.

*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.

^{−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

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

#### 3.2. A tunable directional coupler based on liquid crystals

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]

*l*

_{1}=2

*µ*m,

*l*

_{2}=3

*µ*m,

*l*

_{3}=4

*µ*m,

*l*

_{4}=11

*µ*m. The two waveguides each have an area of 1mm

^{2}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

*n*=1.5 and

_{o}*n*=1.6, respectively. The substrates have a refractive index

_{e}*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

*E*1 and

*E*2 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).

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]

*V*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

*θ*and twist

_{e}*ϕ*angles within the element by using Eq. (29).

_{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

_{e}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]

*xy*-plane. Compared to the finite difference methods reported in [3

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]

*E*] with a vertical linear polarization and wavelength

_{y},E_{z}*λ*=1.5

*µ*m was simulated in the

*z*-direction for a propagation step Δ

*z*=1

*µ*m. The

*E*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

_{y}*E*and

_{y}*E*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.

_{z}*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 7

*V*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.

#### 3.3. Diffraction in a spatial light modulator (SLM)

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

*z*-direction of a plane wave [0,

*E*,0] with vertical linear polarization and wavelength

_{y}*λ*=1

*µ*m through two structures each consisting of 4 pixels as shown in Fig. 11 with dimensions

*l*=4

_{p}*µ*m and

*l*=2

_{p}*µ*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

*n*=1.5 and

_{o}*n*=1.6. For the top left and bottom right pixel, which are in the ON state, the crystalline

_{e}*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).

*E*,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

_{y}*E*component obtained by simulating the propagation of the

_{x}*E*plane wave through the structure over

_{y}*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.

*E*in Fig. 12 show that the original vertical linear polarization of the input plane wave (

_{x}*E*=0) is indeed rotated over 90° after propagation through the ON pixels over 5

_{x}*µ*m to obtain a horizontal linear polarization. As expected, the initial vertical polarization largely remains unchanged in the OFF pixels (

*E*≈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

_{x}*µ*m by 4

*µ*m pixels, the resulting amplitude profile of

*E*consists of several peaks due to elementary Fresnel diffraction of a plane wave. Due to the smaller pixel dimensions, the

_{x}*E*component has a more confined diffraction profile in the 2

_{x}*µ*m by 2

*µ*m pixels. It can be seen from Fig. 12 that due to light diffraction of

*E*in the ON pixels, this component spreads out to partially cover the OFF pixels. This effect is already present for the 4

_{x}*µ*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.

## 4. Conclusions

## Acknowledgments

## 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. |

2. | C. Ma and E. Van Keuren, “A three-dimensional wide-angle BPM for optical waveguide structures,” Opt. Express |

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

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

5. | K. Saitoh and M. Koshiba, “e,” J. Lightwave Technol. |

6. | D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, “Mixed finite element beam propagation method,” J. Lightwave Technol. |

7. | M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. |

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

9. | J. Jin, |

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

12. | G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. |

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 |

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

16. | J. C. Campbell, F. A. Blum, D. W. Shaw, and K. L. Lawlay, “GaAs Electro-optic directional coupler switch,” Appl. Phys. Lett. |

17. | P. G. de Gennes and J. Prost, |

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

19. | E. Buckley, “Holographic Laser Projection Technology,” SID Int. Symp. Digest Tech. Papers |

**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: Year | Journal | Reset

### References

- 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]
- 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]
- 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]
- 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]
- 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]
- D. Schulz, C. Glingener, M. Bludszuweit, and E. Voges, "Mixed finite element beam propagation method," J. Lightwave Technol. 16, 1336-1342 (1998). [CrossRef]
- 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]
- 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]
- J. Jin, The finite element method in electromagnetics, 2nd edition (Wiley, New York US, 2002).
- 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]
- 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]
- G. R. Hadley, "Multistep method for wide-angle beam propagation," Opt. Lett. 17, 1743-1745 (1992). [CrossRef] [PubMed]
- GiD, the personal pre and post processor, http://gid.cimne.upc.es/.
- 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]
- 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]
- 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]
- P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, 2nd edition (Clarendon, Oxford UK, 1993).
- 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]
- 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.

### Figures

Fig. 1. |
Fig. 2. |
Fig. 3. |

Fig. 4. |
Fig. 5. |
Fig. 6. |

Fig. 7. |
Fig. 8. |
Fig. 9. |

Fig. 10. |
Fig. 11. |
Fig. 12. |

« Previous Article | Next Article »

OSA is a member of CrossRef.