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], finite difference (FD) schemes [2, 3, 4] and finite element (FE) [5, 6] implementations. Employing a finite element scheme offers high accuracy and the flexibility to model arbitrarily curved structures. Hybrid edge/nodal finite elements [7] 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, 5] 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:

where k ₀ is the wavenumber in vacuum, ε̿ is the relative permittivity tensor and µ̿ is the relative permeability tensor. Perfectly matched layers (PML) for anisotropicmedia [8] for reflectionless absorption of electromagnetic waves can be incorporated in Eq. (1) by interchanging µ̿⁻¹ and ε̿ with p̿ and q̿ respectively. These are permittivities and permeabilities modified to include the absorption coefficients in the eight PML regions [5] 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.

Download Full Size | PDF

Using a slowly varying envelope approximation, the electric field Ē can be separated into a slowly varying complex field $\overline{\Phi }$ (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 ₀ n ₀ z), where z is the direction of light propagation:

with n ₀ 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:

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) N_x and N_y and quadratic normal (QN) L shape functions for interpolation of the transverse and longitudinal field [7]. 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] is applied to the wave equation from Eq. (1). Substituting Eqs. (2) and (3) into Eq. (1) leads to the beam propagation expression:

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]. Assuming the variation in ∂p_ii/∂z is sufficiently small, the submatrices in Eq. (4) can be evaluated as:

where the coefficients related to the transverse anisotropy in permittivity (q_xy and q_yx) appear in [A_tt] and the coefficients related to the longitudinal anisotropy in permittivity (q_xz, q_zx, q_yz and q_zy) appear in [C_tz] and [C_zt]. The ϕ_z component can be eliminated from the second row of Eq. (4):

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

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:

where A and B are defined as:

with [1] the identity matrix having the same dimensions as the submatrices [A ₁₁], [A ₁₂] and [B ₁₁] which are defined as:

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:

with i the iteration number and

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,6], 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] 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, 6, 12] 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:

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:

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 ₀ and z ₀+Δz to calculate the central derivative of the transverse field with respect to z at z ₀. This derivative can then be inserted into the propagation vector at z ₀ before the wide angle propagation step is carried out towards z ₀+Δ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] 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].

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 d_i=1.5µm (i=1 to 4), the stability parameter ϑ=1 and maximum attenuation α_max=12.5 in the PML regions [5]. 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] 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 E_x and E_y 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 E_z component of the input optical field shown in Fig. 2(d) is calculated from Eq. (12) where the initial derivative of the transverse field $\genfrac{}{}{0.1ex}{}{{\mathrm{d\varphi }}_t}{\mathrm{dz}}\mathrm{at}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 n_x=1.5, n_y=1.6 and n_z=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 E_x and E_y decrease slowly during propagation. Furthermore, E_x and E_y will propagate at different phase velocities c/n_x resp. c/n_y through the medium, which leads to a phase delay Γ=2π/λ(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_z 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) E_x component, (c) E_y component, (d) E_z component.

Download Full Size | PDF

 

Fig. 3. Evolution of E_z upon propagation through the uniaxial layer: (a) d=1.25µm, (b) d=2.5µm, (c) d=3.75µm, (d) d=5mm.

Download Full Size | PDF

For linearly polarized light there is a symmetry axis through the lobes of the E_z 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 E_x and E_y 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 E_z in Fig. 3(a). For d=2.5µm, the E_x and E_y 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, E_x and E_y 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] 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̄ⁱ 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̄ⁱ 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 (Eⁱ_y≠0,Eⁱ_x=Eⁱ_z=0), this yields:

The magnetic field H̄^t of the plane wave transmitted to the anisotropic medium is oriented parallel to the x-axis (H^t_x≠0,H^t_y=H^t_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:

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 n_e=1.6, respectively. The input field has a gaussian intensity profile for E_y, linear vertical polarization (E_x=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:

with ε _⊥=n ² ₀, ε _‖=n ² _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:

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.

Download Full Size | PDF

 

Fig. 5. Intensity profile of the gaussian optical field (a) at z=0µm and (b) at z=10µm.

Download Full Size | PDF

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.

Download Full Size | PDF

As expected, the numerical error on the beam displacement decreases with the number of elements in Fig. 6. This error is about 10⁻⁵ when 300 elements are used and below 10⁻⁶ for 500 elements. The numerical accuracy of the perfectly matched layers that are applied is evaluated in [5]. 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]. 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 ₁=2µm, l ₂=3µm, l ₃=4µm, l ₄=11µm. The two waveguides each have an area of 1mm² and they are separated by a distance d=1µm. To realize a tunable coupler, the waveguides with refractive index n ₁=1.7 are surrounded by a layer of uniaxial liquid crystal with ordinary and extraordinary refractive index n_o=1.5 and n_e=1.6, respectively. The substrates have a refractive index n ₃=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 ε ₁ and ε ₃. 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] 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] 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.

Download Full Size | PDF

 

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.

Download Full Size | PDF

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] which requires all directors to be aligned in the xy-plane. Compared to the finite difference methods reported in [3,4,18] 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,E_y,E_z] with a vertical linear polarization and wavelength λ=1.5µm was simulated in the z-direction for a propagation step Δz=1µm. The E_y 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 E_y and E_z 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 E_y (top) and E_z (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.

Download Full Size | PDF

 

Fig. 10. Evolution of the electric field components E_y (top) and E_z (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.

Download Full Size | PDF

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] 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,E_y,0] with vertical linear polarization and wavelength λ=1µm through two structures each consisting of 4 pixels as shown in Fig. 11 with dimensions l_p=4µm and l_p=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 n_o=1.5 and n_e=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.

Download Full Size | PDF

If we consider the propagation of the plane wave [0,E_y,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 E_x component obtained by simulating the propagation of the E_y 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 E_x in Fig. 12 show that the original vertical linear polarization of the input plane wave (E_x=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 (E_x≈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 E_x consists of several peaks due to elementary Fresnel diffraction of a plane wave. Due to the smaller pixel dimensions, the E_x 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 E_x 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 E_x component obtained after propagation of a plane wave [0,E_y,0] over 5µm through the 4-pixel structure for (a) 4µm by 4µm and (b) 2µm by 2µm pixels.

Download Full Size | PDF

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.