Accurate modelling of the optics of high resolution liquid crystal devices: a reconfigurable liquid crystal grating

ABSTRACT

An accurate method to model the optical behaviour of liquid crystal (LC) devices, particularly suited to devices where diffractive effects are present is described here. An accurate electromagnetic modelling programme that takes into account the full non-uniformity and anisotropy of the LC has been developed. This is combined with an existing in-house LC finite element modelling programme based on the Landau – De Gennes theory, that uses the order tensor representation of the LC orientation and allows an accurate descriptions of structures containing LC defects and small features. The electromagnetic model is based on the total field/scattered field (TF-SF) approach to electromagnetic scattering problems and is implemented using finite differences in the frequency domain (FDFD) in a form that can accommodate perfectly matched layers (PMLs) and periodic boundary conditions. The resultant matrix problem is solved efficiently using an especially adapted form of a sweeping preconditioner and the generalised minimum residual method (GMRes). This method has been implemented in 2D and is demonstrated here with the design and analysis of a reconfigurable blazed phase grating that utilises an LC defect to produce an abrupt fly-back, with the capability of short periods and high diffraction efficiency.

GRAPHICAL ABSTRACT

KEYWORDS:

• Liquid crystal modelling
• diffractive LC device modelling
• LC optics modelling
• liquid crystal gratings

1. Introduction

In comparison with other commonly used optical materials, liquid crystals (LC) have a large birefringence that can be controlled with relatively low electric fields, making them very attractive for many photonic, microwave and terahertz applications. In photonics, apart from traditional displays, they are ubiquitous in diffractive devices, programmable holograms and spatial light modulators, with applications in communications and many other fields [1–7] and more recently, LC geometric phase lenses for beam steering, microscopy, laser tweezing and LiDAR [8–11]. Diffractive LC devices are also the basis for new display systems such as head-up, holographic, diffractive and 3D displays [12], and in general, applications where phase modulation is necessary, something that is not possible with emissive systems. Another range of applications is in shaping and controlling electromagnetic wave propagation in photonics, microwave and terahertz waveguiding devices [13–17].

In many of these applications, changes in the material characteristics (LC orientation, presence of defects) and dimensions of the device features and size and spacing of electrodes, occur at length scales comparable to the wavelength used and accurate modelling becomes essential to design and optimise these devices. This includes not only the balance between biasing electric fields and elastic interactions which will determine the LC orientation and then, the permittivity distribution, but also the accurate modelling of electromagnetic wave propagation through this fully anisotropic and non-uniform material.

We describe here a finite difference in the frequency domain (FDFD) method to model electromagnetic wave propagation, combined with a finite element LC model that yields the permittivity distribution for any state of switching of the LC material. A finite element programme based on the Landau-De Gennes theory [18–20] is used to find the distribution of the LC orientation over the LC cell, from which the high frequency or optical tensor permittivity distribution is calculated. This is then imported into the FDFD electromagnetic solver.

The numerical solution of electromagnetic problems of this kind is normally not simple. Fine discretisation of the problem domain is needed and the problem often presents ill-conditioning, which is usually exacerbated by the use of absorbing conditions that are necessary in many problems of interest. All these circumstances combined will cause low convergence rates, long computation times and large memory requirements. To try to overcome these problems, we use here a procedure that combines an adapted form of a sweeping preconditioner [21], within the total field – scattered field (TF-SF) approach to electromagnetic scattering problems [22–26].

The method has been implemented in two dimensions and is demonstrated here with the design and analysis of a reconfigurable blazed phase grating that uses LC defects to introduce an abrupt fly-back and reduced grating period.

Liquid crystal diffraction gratings and spatial light modulators (SLMs) in general, are common component in numerous optical systems. See for example [27,28]. LC-blazed gratings and SLMs are for many applications realised using arrays of LC cells in transmission or commercial liquid crystal on silicon (LCOS) panels for reflective devices, that consist of square pixels with a small inter-pixel gap. The characteristics of the device are determined by the voltages applied to the electrodes in each pixel. In this way, a blazed grating can be realised arranging the voltages in a stepped manner to achieve the desired-phase response [3,29–32], aiming at a piecewise constant phase distribution on the output plane. However, a number of pixels are required for each period of the grating, leading to small deflection angles. In these conditions, defects are unwanted and the inter-pixel crosstalk, although helping to smooth out the phase in the ramped part of the blazed grating, prevents a sharp fly-back. These effects are increased when the pixel size is reduced, which is necessary to achieve larger deflection angles. Here we study how a LC defect can be utilised to introduce a sharp change in LC orientation with the accompanying change in the phase response in the output plane. A preliminary generic design is presented and routes to optimisation are indicated.

2. Theory

2.1. Liquid crystal modelling

The LC modelling used here is based on the Landau – de Gennes theory and describes the local LC orientation and order by the order tensor Q, from where the director and the scalar order parameters can be extracted [18,19,33,34]. The order tensor, or Q-tensor, is symmetric and traceless and can be specified by five scalar values. This formulation, allowing the LC order to vary, is particularly adequate to study LC defects and other situations where the order varies substantially.

A variational approach, implemented with the finite element method in 3D, is used to find the dynamic evolution of the Q tensor distribution that yields stationary values to the Gibbs free energy functional:

(1) $F(\mathbf{Q},\phi )=\underset{\mathrm{\Omega }}{\int }\left\{f_D(\mathbf{Q})+f_B(\mathbf{Q})-f_E(\mathbf{Q},\phi )\right\}d\mathrm{\Omega }+\underset{\mathrm{\Gamma }}{\int }{f}_S(\mathbf{Q})d\mathrm{\Gamma }$(1)

where f_B is the Landau-de-Gennes bulk thermotropic energy that can be written as a power expansion in the scalar invariants of the Q-tensor about the nematic-isotropic transition and truncated to terms of order four:

(2) $f_B(\mathbf{Q})=\frac{1}{2}A\mathrm{t}\mathrm{r}({\mathbf{Q}}^2)+\frac{1}{3}B\mathrm{t}\mathrm{r}({\mathbf{Q}}^3)+\frac{1}{4}C(\mathrm{t}\mathrm{r}({\mathbf{Q}}^2){)}^2$(2)

where A, B and C are material-dependent parameters.

The elastic distortion energy density f_D can be written as an expansion in the gradients of the Q-tensor as:

(3) $\begin{array}{c}{f}_D(\mathbf{Q})=\frac{1}{2}\{L_1{\mathbf{Q}}_{\alpha \beta ,\gamma }{\mathbf{Q}}_{\alpha \beta ,\gamma }+L_2{\mathbf{Q}}_{\alpha \beta ,\beta }{\mathbf{Q}}_{\alpha \gamma ,\gamma }+L_3{\mathbf{Q}}_{\alpha \beta ,\gamma }{\mathbf{Q}}_{\alpha \gamma ,\beta }\\ +L_4{\sigma }_{\alpha \beta \gamma }{\mathbf{Q}}_{\alpha \delta }{\mathbf{Q}}_{\beta \delta ,\gamma }+L_6{\mathbf{Q}}_{\alpha \beta }{\mathbf{Q}}_{\gamma \lambda ,\alpha }{\mathbf{Q}}_{\gamma \lambda ,\beta }\}\end{array}$(3)

where $L_i$ are elastic coefficients that can be expressed in terms of the material’s Frank elastic constants and ${\sigma }_{\alpha \beta \gamma }$ is the Levi-Civita tensor.

The electric energy density f_E, associated to the applied electric fields is given in general by:

(4) $f_E(\mathbf{Q},\phi )=\frac{1}{2}{\epsilon }_0(\mathrm{\nabla }\phi \cdot \overline{\stackrel{ˉ}{\epsilon }}\cdot \mathrm{\nabla }\phi )$(4)

where φ is the electric potential

The surface term in (1) corresponds to the anchoring energy that represents the attachment between the LC and the bounding surfaces. This can be described by [19,34,35]:

(5) $f_S=aTr({\mathbf{Q}}^2)+W_1({\hat{v}}_1\cdot \mathbf{Q}\cdot {\hat{v}}_1)+W_2({\hat{v}}_2\cdot \mathbf{Q}\cdot {\hat{v}}_2)$(5)

The principal axes of anchoring $(\hat{e},{\hat{v}}_1,{\hat{v}}_2)$ are the easy direction and two mutually orthogonal unit vectors such that $\hat{e}={\hat{v}}_1\times {\hat{v}}_2$. The constants W₁ and W₂ define the azimuthal and polar surface anchoring strength and a determines the resultant easy surface order.

The modelling of the LC hydrodynamics is done following the approach by Sonnet et al. [36,37], which is numerically equivalent to the Qian and Sheng formulation [38], a generalisation of the Eriksen and Leslie theory for the variable order case. With this, and not considering translational flow, the Q-tensor evolution is characterised by:

(6) ${\mu }_1\frac{\mathrm{\partial }{Q}_{ij}}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }f}{\mathrm{\partial }{Q}_{ij}}+{\mathrm{\partial }}_k\frac{\mathrm{\partial }f}{\mathrm{\partial }{Q}_{ij,k}}$(6)

where μ₁ is proportional to the rotational viscosity and f is the total energy density.

Once the Q-tensor distribution is found, the permittivity distribution is then calculated as $\overline{\stackrel{ˉ}{\epsilon }}={\epsilon }_0({\epsilon }_{\mathrm{\perp }}\overline{\stackrel{ˉ}{\mathbf{I}}}+\mathrm{\Delta }\epsilon \hat{n}\otimes \hat{n})$ where $\hat{n}$ is the LC director, ${\epsilon }_{\mathrm{\perp }}$ is the relative permittivity calculated in the direction perpendicular to the long axis of the LC molecules and $\mathrm{\Delta }\epsilon ={\epsilon }_{\parallel }-{\epsilon }_{\mathrm{\perp }}$ is the dielectric anisotropy.

The solution of the coupled equations that describe the LC physics and the electrostatic potential is implemented using the finite elements method [18–20], to find the dynamic evolution of the LC orientation and order distribution. The problem domain is discretized with a mesh of first-order elements. Neglecting translational flow in the LC, the variables on the nodes of the mesh are the five parameters that describe the Q-tensor plus the electric potential. The solution process consists of a time stepping procedure that calculates iteratively the electric potential and the Q-tensor. At each time step, the potential distribution is found first from the current value of the Q-tensor, and then the Q-tensor is calculated in an iterative nonlinear Crank-Nicolson scheme using Newton-Raphson iterations to deal with the nonlinearity. These calculations are repeated within the same time step until consistency between electric potential and Q-tensor distributions is achieved before progressing in the time stepping procedure.

To accelerate the solution process, variable time stepping is implemented using an estimate of the time derivative error to calculate an appropriate step size. Adaptive mesh refinement has also been implemented to even-up the error distribution and increase the efficiency of the calculations [34].

2.2. Electromagnetic modelling

The propagation of electromagnetic waves through the LC cell is modelled using finite differences in the frequency domain (FDFD) [25,26,39,40]. The calculation method starts by formulating Maxwell’s curl equations in the form [41,42]:

(7) $\mathrm{\nabla }\times \vec{E}=k_0{\mu }_r\overrightarrow{\tilde{H}}\mathrm{\nabla }\times \overrightarrow{\tilde{H}}=-j{\eta }_0\vec{J}+k_0{\overline{\stackrel{ˉ}{\epsilon }}}_r\cdot \vec{E}$(7)

where $\overrightarrow{\tilde{H}}=-j{\eta }_0\vec{H}$ is the scaled magnetic field and we assume ${\mu }_r$ is scalar and uniform. Taking the curl of the first equation and substituting the second in the right-hand side leads to the Helmholtz equation for the electric field:

(8) $\mathrm{\nabla }\times \mathrm{\nabla }\times \vec{E}-k_0^2{\mu }_r{\epsilon }_r\cdot \vec{E}=-j\omega {\mu }_0{\mu }_{r}J$(8)

Writing this equation in matrix form and assuming uniformity along y for the 2D case gives:

(9) $\left[\begin{array}{c}-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}0\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }z}\\ 0-{\mathrm{\nabla }}_{xz}^{2}0\\ \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }z}0-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}\end{array}\right]\left[\begin{array}{c}{E}_x\\ \hspace{0.17em}\\ E_y\\ \hspace{0.17em}\\ E_z\end{array}\right]-k_0^2\overline{\stackrel{ˉ}{\epsilon }}\cdot \left[\begin{array}{c}{E}_x\\ \hspace{0.17em}\\ E_y\\ \hspace{0.17em}\\ E_z\end{array}\right]=-j\omega {\mu }_0\left[\begin{array}{c}{J}_x\\ \hspace{0.17em}\\ J_y\\ \hspace{0.17em}\\ J_z\end{array}\right]$(9)

where ${\mathrm{\nabla }}_{xz}^2=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}$.

The domain is discretised using a regular rectangular grid and assuming that the general direction of propagation is z, the nodes are numbered following a frontal arrangement based on layers perpendicular to the z-axis. This is an arbitrary choice, a simple addressing strategy and doesn’t restrict the solutions but is convenient if the propagation is predominantly along the z-direction. The discretised form of the equations for the three components of the electric field from (8) are then assembled, forming a vector of unknowns of the form: $\mathbf{f}=(E_x^{(1)},E_y^{(1)},E_z^{(1)},E_x^{(2)},E_y^{(2)},\dots ,E_z^{(N)})$ for each of the nodes of the grid (denoted by the number in the superscript). This arrangement is more convenient than grouping by field component because it yields a more compact sparsity pattern of the resultant matrix. Similarly, a unified vector of the source terms is defined as: $\mathbf{g}=(J_x^{(1)},J_y^{(1)},J_z^{(1)},J_x^{(2)},J_y^{(2)},\dots ,J_z^{(N)})$.

Equation (9)(9) $\left[\begin{array}{c}-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}0\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }z}\\ 0-{\mathrm{\nabla }}_{xz}^{2}0\\ \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }z}0-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}\end{array}\right]\left[\begin{array}{c}{E}_x\\ \hspace{0.17em}\\ E_y\\ \hspace{0.17em}\\ E_z\end{array}\right]-k_0^2\overline{\stackrel{ˉ}{\epsilon }}\cdot \left[\begin{array}{c}{E}_x\\ \hspace{0.17em}\\ E_y\\ \hspace{0.17em}\\ E_z\end{array}\right]=-j\omega {\mu }_0\left[\begin{array}{c}{J}_x\\ \hspace{0.17em}\\ J_y\\ \hspace{0.17em}\\ J_z\end{array}\right]$(9) becomes now a matrix equation of the form: A f = g.

2.2.1. Mesh truncation

The problem domain and the corresponding mesh need to be terminated with proper boundary conditions. Electric or magnetic walls, where Dirichlet or Neumann boundary conditions apply, as well as periodic boundary conditions, present no implementation problems and bring no undesired consequences to the solution of the resultant matrix problem. If the domain is naturally unbounded, the mesh is usually terminated using an absorbing boundary condition that prevents outgoing waves to be artificially reflected back into the truncated domain. For this purpose, we use the stretched-coordinate formulation of perfectly matched layers (SC-PML) [43–45]. Although the Uniaxial PML formulation [46] is easier to implement, it has been shown [47] that the stretched-coordinate formulation leads to a global matrix with a lower condition number, which in turn leads to a more efficient matrix solution. To apply the SC-PML, the analytic extension $\tilde{x}=x-jf\left(x\right)$ is assumed, where it is customary to choose the frequency-dependent form: $f(x)=\sigma x/\omega$. Setting the coefficient $\sigma$ to zero in the region of interest and with a positive non-zero value in the PML region, the solution of the equations will be an unaltered wave in the region of interest and a wave with exponentially decaying amplitude in the PML. To find the solution, a transformation of coordinates is needed to bring the problem back to the real coordinates: defining $s_x=1-j{\sigma }_x/\omega$, and similarly $s_z$, this is equivalent to transforming the differential operators by:$\widetilde{\frac{\mathrm{\partial }}{\mathrm{\partial }x}}=\frac{1}{s_x}\frac{\mathrm{\partial }}{\mathrm{\partial }x};\widetilde{\frac{\mathrm{\partial }}{\mathrm{\partial }z}}=\frac{1}{s_z}\frac{\mathrm{\partial }}{\mathrm{\partial }z}$.

Applying this transformation to Equation (8)(8) $\mathrm{\nabla }\times \mathrm{\nabla }\times \vec{E}-k_0^2{\mu }_r{\epsilon }_r\cdot \vec{E}=-j\omega {\mu }_0{\mu }_{r}J$(8) leads to a discretised set of equations for the values of the field components at the nodes of the grid given by the matrix problem: A f = g.

The matrix A has a symmetric sparsity pattern and is symmetric Hermitian if the material is isotropic (not the case of LCs) or if it is anisotropic but lossless. The Hermitian property is also lost if PMLs are used.

The TF-SF approach is a convenient method that was initially developed for the FDTD solution of scattering problems [22–26] and consists of defining equivalent sources on an artificial surface surrounding the scatterer in order to replace the external excitation fields.

Explicit use of the vector g in the right-hand side of the matrix equation is avoided using the procedure presented by Rumpf [25,26], especially adapted to FDFD calculations.

This approach is carried out working with the fully assembled matrix and defining a masking matrix R—a diagonal matrix with 1 and 0 elements in the diagonal – to identify the vector components in the scattered field (SF) and the total field (TF) regions. In this form, the source field values that correspond to the scattered field or to the total field regions can be identified easily as: $\mathbf{R}{\mathbf{f}}_{source}$ and $(\mathbf{I}-\mathbf{R}){\mathbf{f}}_{source}$ respectively, where ${\mathbf{f}}_{source}$ is the excitation field. These expressions can then be used to add or subtract the source terms in the equations corresponding to TF and SF regions, resulting in the right-hand side vector: g [26]: $\mathbf{g}=(\mathbf{R}\mathbf{A}-\mathbf{A}\mathbf{R}){\mathbf{f}}_{source}$ avoiding the explicit calculation of equivalent sources. The use of the TF-SF approach is particularly useful to make the PML more effective by removing the applied source field incident on PML boundaries and also for the subsequent post-processing of the results.

2.2.2. Solution of the matrix problem

The most time-consuming part of the solution process is the solution of the resultant matrix problem and in consequence the efficiency of the method chosen for this task is important.

For a problem limited entirely by electric or magnetic walls, or by a periodic boundary condition and in the absence of losses, the resultant matrix is Hermitian and there are several adequate methods to solve the matrix problem. However, the presence of losses or the introduction of PMLs will destroy the symmetry and increase the condition number of the matrix. In those cases, a method like GMRes is convenient but the matrix needs to be preconditioned in order to be solved efficiently. In this work, we choose the sweeping preconditioner proposed by Engquist and Ying [21] specially developed to solve matrix problems originated by the numerical solution of the Helmholtz equation. The procedure consists on finding an approximate block LDL^T factorisation of the original matrix A and use this as a preconditioner. However, this is done in a block Gauss elimination process that requires the inverse of the diagonal blocks. The basic idea of the method consists of using a frontal ordering approach to the generation of the matrix A, where each of the blocks correspond to layers in the problem domain. With this, the solution for the Schur’s complement in the block Gauss elimination can be seen as the approximate solution of the problem corresponding to just a number of layers which are bounded by absorbing boundary conditions. Sweeping through the domain a number of layers at the time in this manner, will give an approximation of the desired blocks and allow to construct the preconditioning matrix. Although valid for full 3D problems, this preconditioning procedure has been implemented here for 2D problems in conjunction with the TF-SF approach to scattering problems and is demonstrated next with the analysis of a LC diffraction grating.

3. Example of application: a liquid crystal reconfigurable blazed phase grating

The methods described above are used here to demonstrate a reconfigurable LC blazed phase grating that uses the generation of defects to obtain a sharp fly-back.

Traditionally, LC defects were seen as unwanted features that deteriorate the device performance. However, in many devices their presence and evolution are essential to their operation, as for example in the case of bistable displays [48–53].

We present here a basic design for a reconfigurable LC blazed phase grating where a defect is created and used to introduce an abrupt fly-back. This can reduce considerably the grating period and in consequence, increase the possible deflection angles.

3.1. Transmission case

The first case studied here consists of a cell 10 μm wide and with a cell gap (height) of 3.3 μm. Two small electrodes 0.5 μm wide and separated by 0.5 μm are placed in the centre of one surface while the other is covered completely by a ground electrode. The LC material parameters are chosen for this example as those of 5CB, with the values: K₁₁ = 6.4 pN, K₂₂ = 3 pN, K₃₃ = 10 pN, ${\epsilon }_{\parallel }=13$, ${\epsilon }_{\mathrm{\perp }}=6.7$, n_o = 1.544 and n_e = 1.736 [33,54]. While in conventional LC gratings, the alignment is perpendicular to the grating in order to avoid the disclination lines, in this case it is chosen along the same direction of the grating. Periodic boundary conditions are used at the sides of the cell and a pretilt angle of 5° is applied to the top and bottom surfaces.

When the voltage is applied, pincement occurs in the vicinity of the inter-electrode gap and two disclination lines are created [55]. The −1/2 defect line then moves towards the centre of the cell reaching a stationary position in approximately 20 ms while the +1/2 defect remains close to the corner of the driving electrode.

Depending on the applied voltage, the position of the −1/2 defect and the corresponding field distribution over the output plane vary. The results shown in Figure 1(a) correspond to 5 V applied to one of the bottom electrodes while the other is connected to ground. The top electrode is also fixed to ground. For clarity of illustration, the calculated directors have been here interpolated on a coarse regular grid. Different voltages were tried and it was found that 5 V produced the results that best approximate an ideal blazed grating shape for the phase distribution over the top surface. This also corresponds to the case were the −1/2 defect migrates to the position closer to the middle of the cell gap.

Figure 1. (Colour online) (a) Director distribution and electric potential contour lines when 5 V are applied to one of the bottom electrodes. the circle shows the stationary position of the −1/2 defect, the +1/2 defect remains stationary just above the bottom surface. (b) Field distribution: magnitude of E_x component when a plane wave is normally incident on the bottom surface of the LC layer. Dimensions are in microns; x-z cross-section is shown.

From the director distribution results, the permittivity tensor can be calculated and interpolated onto a regular grid suitable for the FDFD calculation of the electromagnetic response. For this purpose too, the structure is extended at the top and bottom adding a narrow air gap and a perfectly matched layer, with periodic boundary conditions imposed on the sides. A plane wave excitation of wavelength 550 nm and polarised along the x-direction is applied at the bottom of the LC layer (z = 0 in Figure 1(b)). The combination of wavelength and cell gap was chosen to obtain approximately a 2π phase range in the grating. Figure 1(b) shows the magnitude of the x-component of the electric field over the whole cross-section. There is a slight polarisation conversion, predominantly in the defect region but the y-component of the electric field remains low over the top surface.

Figure 2(a) shows the resultant-phase distribution in the output plane, extracted just above the top surface of the LC layer. The far-field diffraction pattern of a single cell is shown in (b) and the diffraction pattern of an array of 100 cells in (c). These results indicate that just with two electrodes the basic blazed grating shape can be approximated, and suggest that a more versatile structure consisting of a continuous array of small, regularly spaced electrodes on the bottom surface, which could then be connected to appropriate voltages, can be used to improve the approximation or to change the period of the grating, or in general, to approximate other desired phase distributions on the top surface or output plane.

Figure 2. (Colour online) (a) Phase distribution on the top surface when a plane wave is incident normally on the bottom surface, (b) Far-field diffraction pattern of a single cell (solid blue), (c) Far-field diffraction pattern of an array of 100 cells. the dotted red lines indicate the shapes corresponding to an ideal blazed grating of the same period and total phase range.

3.1.1. Regularly spaced electrodes – ten electrode case

We consider now a section of 10 µm with 10 regularly spaced electrodes on the bottom surface of the same dimensions and spacing as those in Figure 1. The voltages applied to the central electrodes are 0 and 5 V as those in Figure 1 and the remaining electrodes were initially set to the voltages corresponding to the potential values at the points corresponding to the centre of their position in the two-electrode structure, as the start of an iterative optimisation procedure.

A full automatic optimisation of the voltages, involving repeated use of the LC and electromagnetic modelling could be prohibitively onerous if too many iterations are needed. However, a simplified method as that used in [31] can make this viable.

For this example, only a simple manual iterative procedure of successive corrections was performed, consisting of changing slightly the voltages at the electrodes in order to improve the shape of the phase distribution. Satisfactory results are achieved after only four or five iterations.

The director distribution and the magnitude of the E_x component are shown in Figure 3. The two central electrodes are set at zero and 5 V as in the previous case, while the others are set as indicated in Figure 4(a) after a simple optimisation. Figure 4(b) shows the phase distribution in the output plane (top surface of the LC layer). The corresponding far-field diffraction pattern for an array of 100 cells appears in Figure 4(c).

Figure 3. (A) (Colour online) Director distribution for the 10-electrode device. (b) Magnitude of the E_x component of the field when an x-polarised plane wave is incident normally from the bottom. the LC region is the central region shown above and is topped by small air layers and PMLs at both ends.

Figure 4. (Colour online) (a) Voltage pattern and values on the electrodes, (b) E_x component phase distribution on the output plane. the red dashed line indicates the shape corresponding to an ideal blazed grating of the same period and phase range, (c) Far-field diffraction pattern of an array of 100 cells.

For the array of 100 cells, the diffraction efficiency for the −1 order mode is 67.8% and the ratio between the 0-order and the −1 orders of diffraction is 31.8%.

The presence of the defect causes a certain amount of twist to occur in the central region of the cell and this causes a degree of polarisation conversion, which is illustrated in Figure 5(a). The solid line corresponds to the magnitude of the E_x component of the field when an x-polarised plane wave is used as input and the red dashed line shows the magnitude of E_y, the leaked field. This effect can be reduced substantially when a LC material with a higher value of the twist elastic constant, K₂₂, is used. Figure 5(b,c) show the field component values at the output plane when the twist elastic constant is increased by a factor of 1.2 and 1.5, respectively.

Figure 5. (Colour online) (a) Amplitude distribution of the electric field components over the output plane for 5CB, E_x solid line, E_y dashed line. (b) Amplitudes if K₂₂ is increased by a factor of 1.2. (c) Amplitudes if K₂₂ is increased by a factor of 1.5.

Changing the twist elastic constant in this form, the polarisation leakage is drastically reduced but there is almost no appreciable effect on the amplitude or phase distribution of the E_x component of the field or the resultant diffraction pattern.

As seen in Figure 5, the amplitude of the E_x component of the field varies significantly through the cell. However, there is little effect of this variation on the far field diffraction pattern. Only a small improvement is observed if the diffraction pattern is calculated assuming a uniform amplitude distribution and this difference is because, for simplicity, the basic optimisation procedure used was looking for the best approximation to the phase response rather than the resultant diffraction pattern.

Apart from seeking a material with higher K₂₂ to improve the performance, a higher dielectric anisotropy would be beneficial since this would lead to a thinner cell and consequently, more tolerance to variations on the incidence angle and more uniformity in the output electric field amplitude.

3.1.2. Changing the grating period

To illustrate the reconfigurabilty of this grating, the period of the structure is now reduced from 10 µm to 8 µm, comprising eight electrodes and keeping the dimensions and spacing of the electrodes unchanged. The two central electrodes are fixed at 0 and 5 V and the rest of the voltages, shown in Figure 6(a), are found by the simple procedure of successive corrections described earlier. Again in this case, a defect pair nucleates and the −1/2 defect migrates to the middle of the cell allowing the overall abrupt change in the director orientations at either side.

Figure 6. (Colour online) (a) Voltage values on the electrodes. (b) Phase distribution of the E_x component of the field on the top surface. (c) Far-field diffraction pattern of an array of 100 cells. Dashed line indicates the shape corresponding to an ideal blazed grating of the same period and phase range.

The results are shown in Figure 6. The deflection angle of the −1 diffraction order increases from 3.16° to 3.96°, with a diffraction efficiency for the −1 order mode of 63.5% and a 0-order to −1 order ratio of 41% for an array of 100 cells.

The phase distribution over the output plane, Figure 6(b), and also the diffraction pattern, Figure 6(c), depart from the ideal a bit more than what is observed in the ten electrode case, but this is expected since less elements are used to approximate the ideal response.

3.2. Reflection case

A common and convenient form of implementation of devices like those described above is using LCOS technology. In this case, the device operates in reflection, the excitation is from above and the bottom electrodes are directly connected to the CMOS layer. As a difference to the standard LCOS devices, where the pixel electrodes are themselves the reflector, in this case, a separate reflection layer would be required.

Considering this configuration, the structure chosen is 10 µm wide by 1.7 µm tall and is filled with 5CB, the same LC material as in the structure used in transmission. The cell gap has been adjusted to have approximately $2\pi$ phase range.

The top surface of the LC layer is covered by a continuous ground electrode and the bottom surface has an array of electrodes of the same dimensions and spacing as before. One electrode, marked with the red outline in Figure 7(a), is connected to 3 V, other (solid black) is connected to ground and the rest, shown in grey, are set to voltages found with the iterative procedure described before, starting from the potential distribution originated by the two central electrodes. Again in this case, just a few iterations are needed to improve considerably the phase distribution and the corresponding diffraction pattern.

Figure 7. (Colour online) (a) Director distribution and equipotential lines over the LC cell when 3 V are applied to one of the bottom central electrodes (marked with red outline) and other is grounded (solid black); the rest of the voltages are chosen to approximate the ideal phase distribution (grey). (b) Magnitude of E_x component when a downwards propagating plane wave is normally incident on the top surface of the LC layer and reflected at z = 0. Dimensions are in microns; xz-plane shown.

The resultant stationary director distribution is shown in Figure 7(a), where the −1/2 defect has already reached the final position near the top of the cell. To study the light propagation through the cell, the model is extended adding an air region at the top and terminating the structure above with a PML. At the bottom of the LC layer, a reflective layer is modelled as a dielectric with large negative imaginary permittivity (a value of 10,000 was used here). Figure 7(b) shows the distribution of the E_x component of the field over the cross-section when a plane wave propagating downwards is incident at z = 1.7 µm. An air layer and a PML were added at the bottom but this is not really necessary since the reflector at z = 0 practically stops all transmission through. The output-phase distribution, shown in Figure 8(a), is extracted at a position in the top-air gap, just above the excitation plane at z = 1.7 µm.

Figure 8. (Colour online) (a) Voltages on the electrodes. (b) Phase distribution of the E_x component of the field on the top surface. the red line indicates the phase corresponding to an ideal blazed grating of the same period and phase range. (c) Far-field diffraction pattern of an array of 100 cells.

The phase distribution over the output plane in the air region, slightly above z = 1.7 µm and the corresponding diffraction pattern are shown in Figure 8(b,c). The resultant diffraction efficiency for the −1 order mode is 76.4% and the 0-order to −1 order ratio is 6.6%. The amplitude of the output field varies in a similar way as in the transmission case using the parameters for 5CB, but with a slightly reduced total variation as shown in Figure 9. However, the thinner cell does not permit twist to occur as easily as in the transmission case and the leakage to the y-polarisation is negligible, even with the low value of K₂₂ of 5CB.

Figure 9. (Colour online) Amplitude distribution of the electric field components over the output plane for 5CB in reflection cell, E_x solid line, E_y dashed line.

The amplitude variation over the output plane is no impediment to achieve a good far field response as shown by the diffraction pattern in Figure 8(c). This is achieved here with a grating period of 10 µm with 10 electrodes per period. As with conventional LCOS gratings, less electrodes per period are possible, but with a reduction in the fidelity of the phase response, as shown earlier.

Conventional LCOS SLM pixels are limited in size due to the strong fringing field effects at small pitch. With a typical small size pixel of about 4 µm, a grating with 10 levels would have a period of 40 µm while with the proposed structure the period would only be 10 µm long, resulting in larger deflection angles. At the same time, while in a conventional LCOS device, the fly-back region would be of the order of 4 µm to 5 µm for the same phase range [32], the fly-back region of the structure described above is only of 1.1 µm.

4. Conclusions

A procedure to model accurately the optical behaviour of LC cells has been presented. It combines in this implementation, an accurate LC modelling method, which is particularly suitable to study structures containing defects and small features, with an efficient electromagnetic solver based on a FDFD formulation of the electromagnetic problem. This leads to a matrix problem that is solved using GMRes. Electromagnetic problems of this type are usually ill-conditioned and this is here treated for a 2D implementation with a modified form of a sweeping preconditioner.

The electromagnetic solver requires the values of the dielectric permittivity in each of the nodes of the FD grid and these are supplied by the LC model.

Using the same strategy followed by this procedure, a less computer-demanding constant-order LC solver, based on the Oseen-Frank formulation, can be substituted here to study problems where order variation can be neglected, as, for example, in microwave and terahertz LC devices, where diffraction effects are important but LC order variation is not.

To illustrate the capability of the method, a 2D case of a reconfigurable LC grating is proposed and analysed. It is shown that a structure consisting of a number of regularly spaced small electrodes can be configured into a periodic array where the electrode voltages are selected to shape a desired phase distribution on the output plane. In this case, this is shown for a blazed grating where a LC defect is utilised to produce the abrupt phase profile. Good results are shown for 5CB, a common LC material but not particularly suitable for this application. It is also shown that a careful choice of material can lead to thinner cells and lower polarisation leakage. Optimisation of the structure is also possible and even with a simple procedure results can be improved substantially in a few iterations.

Built in this form and using LC defects to produce the abrupt fly-back, the period of the grating can be drastically reduced compared with the standard implementation of LC gratings with LCOS SLMs, designed to approximate uniform levels of amplitude and phase over each pixel.

Disclosure statement

No potential conflict of interest was reported by the author(s).