The exact molecular reorientation model for nematic liquid crystals taking into account all diagonal Frank elastic constants and using two angles to describe director orientation is presented. Solutions and simplified equations are shown for the most common planar and chiral configurations. Gaussian beam propagation simulated using fully vectorial Beam Propagation Method in nonlinear case is also provided. Detailed comparison between exact solutions and single Frank constant approximation is made. However, no significant differences between these two models were found neither in beam propagation nor in polarization distribution, some difficulties may occur in choosing single Frank constant especially when it comes to quantitative results. Presented results correspond to a propagation of a beam of the Gaussian or topologically similar shapes.
© 2012 Optical Society of America
Nematic liquid crystals are unique materials with many applications, like stress and temperature sensors, tunable filters, lasers and the most important and widely used - displays [1–3] where the nematic cell is controlled by the external electric field. Liquid crystals may be also controlled by the optical signals due to the nonlinear properties of these materials. Light beam launched into the cell can become non-diffractive and form a nematicon (i.e. a soliton in nematic liquid crystals) [4–6]. This phenomenon was observed in many experiments [7, 8] as well modeled numerically [9–12]. To simulate light propagation in a nematic cell it is needed to model electromagnetic field propagation as well as reorientation of liquid crystals molecules, because in nematics, nonlinearity mainly comes from the reorientation of the molecules. One of the methods is to use Beam Propagation Method (BPM) [13,14] for field propagation and Frank-Osseen molecular model to describe director orientation (i.e. mean direction of molecules). The most common approach to date is to use a single Frank constant approximation and one [15, 16] or rarer two angles [17–19] to describe director orientation. When it comes to chiral nematics (i.e. cholesterics) there is an additional need to use vectorial BPM .
In this paper we present exact equations describing molecular reorientation in chiral and non-chiral nematic liquid crystals based on the Frank-Osseen theory. All Frank elastic constants are kept and two angles to describe director orientation are used. Comparison between this model and one using a single constant approximation is shown. At last we present results of beam propagation obtained with vectorial BPM combined with molecular reorientation model.
At first we deal with molecular reorientation in chiral and non-chiral nematic liquid crystals. The setup and coordinate system are shown in Fig. 1. Two angles to describe director orientation are used so n⃗ = [cosθ, sinθ sinϕ, sinθ cosϕ]. To derive equations describing molecular reorientation we start up with Frank-Oseen equation on free energy density [21–23]:Eq. (1) and then Euler-Lagrange equations of the following form are used: Equations (2) give two following equations describing molecular reorientation: Equations (3) and (5) can be simplified by neglecting derivatives with respect to z. Such assumption means that intensity and polarization slowly varies along the direction of propagation i.e. we use slowly varying envelope approximation. The following equations are obtained:
Another significant simplification can be made by using single Frank constant approximation. We preserve all derivatives with respect to x, y and z directions and assume K11 = K22 = K33 = K. In this case Eqs. (3) and (5) simplify to the already known form [17–19]:Eq. (7) simplifies to the following form: Equation (8) is fulfilled for . Moreover, for planar configuration (G = 0) Eq. (8) is fulfilled for or for . It means that Eq. (11) can only be used for homogeneous field along the y direction or (in planar case) along the x or y direction. In a single constant approximation Eq. (8) vanishes so Eq. (11) can be used without additional assumptions.
In the following part we will discuss only planar configuration (G = 0). When we assume ϕ = 0° and Ey = 0 so the reorientation occurs in xz plane Eq. (8) simplifies to the form:Eq. (7) if fulfilled for or for i.e. for homogeneous fields along either x or y direction.
As it was stated, the above Eqs. (11) and (12) can be used for homogeneous fields only. For finite beams reorientation will occur not only in a single plane and will lead to both ϕ and θ reorientations. Such molecular reorientation occurs mainly at the edges of the beam, where the field evanesce and it gives the lowest energy of orientation. In other words reorientation of molecules in a single plane will not lead to the lowest energy state and thus ϕ and θ reorientation have to occur.
Assuming that ϕ = 90° and Ez = 0 i.e. the reorientation occurs in xy plane, Eq. (8) simplifies to:Eq. (7) vanishes for G = 0 the above equation needs no additional assumptions.
2. Numerical results
In this section numerical simulations of molecular reorientation and beam propagation in nematic liquid crystals are presented. For simulations Gaussian beam containing the electric (E⃗) and magnetic (H⃗) field is used:
For modeling molecular reorientation only electric field components of Gaussian beam are used. To obtain desired power (P) of the beam the amplitudes E0x and E0y were chosen to fulfill ∬Szdxdy = P were Sz - is a z component of Poynting vector.
2.1. Molecular reorientation
Calculations of molecular reorientation are made for a cell of L = 25μm in height and W = 50μm in width on a 2D grid of 50x100. Electric field amplitude corresponds to a Gaussian beam of a power P = 100mW. Optical field waist is equal ω0 = 3μm (FWHM≈ 7μm). Parameters of liquid crystals correspond to 6CHBT liquid crystals [24, 25] with anisotropy Δε = 0.47 and Frank elastic constants K11 = 8.8pN, K22 = 3.7pN and K33 = 9.5pN. Equations (7) and (8) are solved using SOR (Succesive Over-Relaxation) iterative method [26, 27], with relaxation parameter Ω = 1.8 and resolution Δx = Δy = Δh = 0.5μm if not otherwise stated. Typically convergence is obtained after no more than 300 iterations which is discussed later. Initial and boundary conditions depend on liquid crystal configuration. We consider chiral and non-chiral (planar) configurations. In planar configuration initial conditions are as follows and ϕ(x,y) = ϕ0. In many simulations reorientation occurs with Freedericksz threshold  therefore molecules (except boundary conditions) where realigned in relation to the initial conditions as follows Δϕ = +0.01° and Δθ = −0.01°. For instance for Ex ≠ 0 and Ey = 0, Ez = 0, terms f(E⃗) and ξ(E⃗) are equal 0 (vanish) so no reorientation can occur. Realigning the molecules makes terms f (E⃗) and ξ(E⃗) non-zero and calculations with a threshold are possible. In chiral configuration initial conditions are ϕ(x,y) = ϕ0 +G′x where and . In calculations pitch of a cholesterics p and pitch induced by the boundary conditions p′ are equal to the cell height so p = p′ = L and G′ = G. For each configuration (chiral and planar) at the edges of the system initial conditions are fixed and constitute the boundary conditions. In the next paragraphs we will discuss solutions of Eqs. (7) and (8) and compare them to the solutions of Eqs. (9) and (10) assuming single Frank constant approximation. We will also present empirical equations which can help to determine value of an elastic constant K when using single constant approximation.
At first let us consider planar configuration where molecules initially lie in the yz plane (θ(x,y) = const. = 90°) and are rotated through some angle ϕ0 corresponding to the z axis. Solutions for ϕ0 = 0° and electric field applied along the y-axis (α = 0°) are presented in Fig. 2. Assuming reorientation only in yz plane Eq. (7) simplifies to Eq. (11) but, as it was stated before, it can be only used for homogeneous field along either x or y direction, because Eq. (8) is satisfied only for homogeneous fields. In fact there is some θ reorientation up to ±2°. Such reorientation can be neglected so the simplified Eq. (11) can be used. Moreover in a single constant approximation there is no θ reorientation (it is below 10−5 deg. which is caused by numerical error). It is consistent with our previous conclusions made while describing Eq. (11). However simplified equation can be used it is difficult to determine a proper value of K in a single constant approximation. In the considered case reorientation occurs only in yz plane. Due to our simulations performed for different liquid crystals and setups empirical equation can be used to determine elastic constant:(16) gives well reproduction of ϕ reorientation not only in discussed setup but in all the following configurations including chiral one. However, maximum values of ϕ reorientation are comparable between both sets of equations, analyzing results obtained with exact Eqs. (7) and (8) shows that the shape of ϕ reorientation is elliptical in contrast to the results obtained using Eqs. (9) and (10) where the reorientation is rather circular. Elliptical shape is caused by the differences between Frank constants corresponding to different deformations (i.e. K11 - splay, K22 - twist and K33 - bend). It was also checked that the elliptical shape of reorientation is neither caused by the boundary conditions nor by using two angles to describe director orientation. Even in one-angle model while using all Frank elastic constants the reorientation has elliptical shape.
When assuming initial conditions ϕ0 = 0° and θ(x,y) = 90° as in the previous case and an electric field applied along the x-axis (x-polarized α = 90°), reorientation will mainly occur in the xz plane (ϕ = const.). Solutions of this configuration are presented in Fig. 3. In this case Eq. (8) simplifies to Eq. (12). However, Eq. (7) has solutions only for homogeneous fields along either x or y directions, it appears that ϕ reorientation is relatively weak (up to ±5°), corresponding to θ reorientation (more than 70°) so the simplified Eq. (12) may be used. However, maximum θ reorientation obtained using Eqs. (7)(8) and Eqs. (9)(10) is the same, the shape of reorientation in a single constant approximation is elliptical in contrast to the circular shape obtained with full equations. To determine K in this configuration we propose an equation:
Consider another planar configuration where the molecules lie in yz plane (θ(x,y) = 90°) and are aligned along the y-axis (ϕ0 = 90°) and electric field is applied along the x-axis (i.e. is x-polarized α = 90°). Numerical results are shown in Fig.4. The reorientation is expected to occur in xy plane only so Eq. (8) can be simplified to Eq. (13). For planar configuration Eq. (7) vanishes so Eq. (13) is an exact one. In fact, results given by exact Eqs. (7) and (8) as well as simplified Eqs. (9) and (10) show that there is no reorientation in yz plane (ϕ reorientation is below 10−5 deg.). To determine value of K for this configuration we propose the following equation:Fig. 4 it is visible that not only maximum values of reorientation but also a shape is the same for both (exact and simplified) sets of equations.
In the previously analyzed, particular planar configurations, equations could be simplified and differences between results obtained using single constant approximation and full equations were not so significant. Moreover, it was quite easy to determine proper K value as the reorientation occurred in a single plane. For now, let us consider more general case, where molecules initially lie in the yz plane and are rotated through an angle ϕ0 = 45° and the applied field is polarized at α = 45° (Ex = Ey). In such case both ϕ and θ reorientations are important and need be taken into account. Solutions are presented in Fig. 5. However, maximum values of reorientation obtained with full Eqs. (7) and (8) and equations using single constant approximation (9) and (10) are comparable, the shape of ϕ reorientation is not only elliptical but also the long axis of ellipse does not lie along the y direction. This rotation of the long axis vanishes in one-angle model, while using only angle ϕ to describe director orientation.
In this configuration, where as well ϕ as θ reorientation is important, it is very difficult to chose proper K value.
Another simulations were made for chiral configuration where molecules are rotated through an angle of 360°. Optical field is polarized at α = 45°. Solutions are presented in Fig. 6. Results obtained with full equations and simplified equations (using single constant approximation) are nearly the same. Total reorientation and the shape of reorientation are comparable. To determine proper K value we propose to use Eq. (16).
Convergence, timing and numerical uncertainties
To determine convergence of Eqs. (7) and (8) we have performed simulations of planar configuration where ϕ0 = 45°, w0 = 3μm, α = 45°, P = 100mW, Ω = 1.8, K11 = 8.8pN, K22 = 3.7pN, K33 = 9.5pN. Calculations were performed as previously for a cell of L = 25μm in height and W = 50μm in width. We have measure maximum positive ϕ reorientation and maximum negative θ reorientation as only these occur in such configuration. Reorientation versus number of iterations is presented in Fig. 7. Convergence of the analyzed equations is very well, but what has to be stressed it strongly depends on resolution. For instance for Δh = 0.5μm (grid of 100x50 points) convergence is obtain after approximately 50 iterations, but for Δh = 0.1μm convergence is not obtained even for 300 iterations. We have also compared convergence of Eqs. (7) and (8) with Eqs. (9) and (10) (see Fig. 7). Convergence of both sets of equations is comparable.
Calculations were made on PC with an Intel Pentium IV clocked at 2.66 GHz and 1.5GB of RAM. Computation time of the above planar configuration on a grid of 50x100 points with 300 iterations was up to 7 sec. for full equations and 6.5 sec. for simplified ones. For a grid of 100x200 points timing was 28 sec. and 26.5 sec. respectively.
To estimate the numerical uncertainty calculations were made for a planar configuration with no external field applied and a noise of amplitude Anoise = 45° added to the initial conditions. Uncertainty was estimated to be not grater than 0.3°.
To show the changes in the results due to different numerical steps, calculations for planar configuration were performed. Results of 3000 iterations performed on a grid of 50x100 (Δh = 0.5μm), 100x200 (Δh = 0.25μm) and 250x500 (Δh = 0.1μm) points using both models were compared. Differences in maximum and minimum reorientation between models did not exceed 2°. Differences were mainly caused by choosing slightly inappropriate K value. When it comes to results for different resolution there was no change between the results (with accuracy ± 0.1°).
2.2. Beam propagation
In this section we describe light propagation in nematic liquid crystal cell. To model molecular reorientation Eqs. (7) and (8) have been used. For modeling optical field propagation we use Full Vector Beam Propagation Method (FVBPM) [9, 20] derived directly from Maxwell equations with no paraxial approximation. Nevertheless, paraxial approximation is indirectly preset as we neglected derivatives with respect to z in molecular reorientation model. The following equations were used to describe light propagation:
Equations (19) were rewritten using finite differences and solved using Runge-Kutta 4th order method. Calculations were performed on a 2D transverse grid with resolution Δx = Δy = 0.5μm. The same resolution was taken to calculate molecular reorientation. Longitudinal resolution was set to Δz = 0.01μm. Molecular reorientation was recalculated after each 100nm of propagation. At the boundaries Dirichlet conditions were used. Simulations were performed for two different liquid crystals 6CHBT (K11 = 8.8pN K22 = 3.7pN K33 = 9.5pN, Δε = 0.44@793nm) and 1110 (K11 = 21pN K22 = 8pN K33 = 21pN Δε = 0.13@1064nm). Gaussian beam described by (14) was launched into the system. Light was propagated along the z axis. Simulations showing beam propagation in a cell filled with liquid crystals in a planar configuration are shown in Fig. 9, 10 and 11. The beam is polarized at α = 45°. Comparison between results obtained using full Eqs. (7) and (8) and results obtained using single constant approximation (Eqs. (9) and (10)) is also made. There are no significant differences neither in beam propagation nor in its polarization distribution between these two models. However, it has to be stressed that it is very important to choose appropriate K value in a single constant approximation. In Fig. 9 and 10 it appears that for inappropriate value of K = 19pN (right) light distribution and beam polarization differs from the exact solution (left). In the xz plane the beam ”breathes” stronger in a single constant approximation. It is caused by the fact that in a single constant approximation ϕ reorientation along the x and y direction is comparable and nearly symmetrical in contrast to the non-symmetrical reorientation appearing when using exact equations (see Fig. 5). Moreover for inappropriate K value the portion of the x-polarized light propagated nearly on axis is larger, and the x-polarized light guided by the y-polarized soliton propagates on a greater distance (see. Fig. 9). It has to be stressed that in a single constant approximation changing K value leads to the same results as changing power. So in fact, choosing appropriate K value in a single constant approximation is about finding the right power values. In Fig. 9 and 10 at the end of propagation the extraordinary part of light slightly changes direction of propagation. It is caused by the boundary which is very close to the propagated beam and has fixed anchoring conditions. Simulations performed for wider numerical grids (in y direction) show no such effect. To compare exact equations with single constant approximation walk-off angle versus input power was also examined (Fig. 8). Both models give comparable results with differences not greater than 0.3°.
When it comes to higher-birefringence liquid crystals 6CHBT, the results obtained with both models are comparable when using correct K value. However, there are some differences comparing to the low birefringence liquid crystals 1110. Obviously, the minimum power needed to form a soliton is significantly lower (about P ≈ 1mW) but there is also a change in polarization propagation. For 1110 some part of the x-polarized is propagated nearly on axis and the rest is guided by the y-polarized soliton. In 6CHBT whole x-polarized light is propagated straight ahead and only the y-polarized forms a soliton.
Simulations of light propagation in a cholesteric configuration are shown in Fig. 12. Cholesterics is rotated through an angle of 360°. The beam is y-polarized. Equations and results are much more stable for y-polarized light in comparison to α = 45°. It may be connected with paraxial approximation which is better fulfilled for the y-polarized light because there is mostly ϕ reorientation and there are nearly no changes in polarization along the z axis. It is visible that for appropriately chosen K the result obtained with both models are comparable. However, diffraction of the beam is confined in yz plane, in the xz plane some ”breathing” appears. It is caused by the presence of the layers in cholesterics. Light reflects from the edge of the layer and since there is no distinguished direction, intensity distribution is symmetric with respect to the center of the cell.
Summarizing, we have presented exact equations describing molecular reorientation in nematic liquid crystals. We have taken into account all derivatives with respect to all directions x, y and z and all diagonal Frank constants K11,K22, K33. To the best of our knowledge it is the first time these equations have been explicitly written. We have shown similarities and differences between molecular reorientation obtained using full equations and equations with single Frank constant approximation. Most of the results obtained using these two models are comparable, but it is crucial to choose an appropriate K value. We also proposed a simple, empirical equation which facilitates calculating K in some simple cases. Full equations neglecting only derivative with respect to z were combined with FV-BPM to provide beam propagation in a liquid crystal cell. Calculations were performed for planar and chiral configurations. No significant differences were observed between results obtained using all Frank constants in comparison to single constant approximation. However, in many cases it is very difficult to choose a correct K value. In such cases we propose to recalculate molecular reorientation using full equations and then recalculate reorientation for a few different K values (using single constant approximation) typically fulfilling K22 < K < K33 and found one which gives results close to the results obtained using full equations. Such K value can be used to calculate molecular reorientation and be combined with FV-BPM. Such approach will significantly shorten computation time as the Eqs. (9) and (10) are much simpler than full equations (7) and (8). Since qualitative results obtained using single constant approximation are similar to the results obtained using all Frank elastic constants, choosing appropriate K value is a matter of scaling the power at which particular phenomena occurs. In other words changing K is the same as changing power P so for qualitative results it is not so crucial and nearly any K value can be assumed. Finally it has to be stressed that presented results and conclusions correspond to the beams of a Gaussian or topologically similar shapes. For other cases, for instance vortexes or other topological defects, results may be different but this is out of scope of our analysis.
This work was financially supported by the National Science Centre. Filip Sala was supported by a scholarship from the Center for Advanced Studies.
References and links
1. D. W. Berreman and W. R. Heffner, “New bistable cholesteric liquid-crystal display,” Appl. Phys. Lett. 37, 109–111 (1980). [CrossRef]
2. S.-Y. Lu and L.-C. Chien, “A polymer-stabilized single-layer color cholesteric liquid crystal display with anisotropic reflection,” Appl. Phys. Lett. 91, 131119 (2007). [CrossRef]
3. B. Bahadur, Liquid Crystals: Applications and Uses (World Scientific Publishing Co. Pte. Ltd., 1995).
4. G. Assanto, M. Peccianti, and C. Conti, “Nematicons: Optical spatial solitons in nematic liquid crystals,” Opt. Photonics News 14, 44–48 (2003). [CrossRef]
5. A. Alberucci, A. Piccardi, M. Peccianti, M. Kaczmarek, and G. Assanto, “Propagation of spatial optical solitons in a dielectric with adjustable nonlinearity,” Phys. Rev. A 82, 023806 (2010). [CrossRef]
6. G. Assanto and M. A. Karpierz, “Nematicons: self-localised beams in nematic liquid crystals,” Liq. Cryst. 36, 1161–1172 (2009). [CrossRef]
9. F. A. Sala and M. A. Karpierz, “Numerical simulation of beam propagation in a layer filled with chiral nematic liquid crystals,” Photon. Lett. Pol. 1, 163–165 (2009).
10. A. Alberucci and G. Assanto, “Propagation of optical spatial solitons in finite-size media: interplay between nonlocality and boundary conditions,” J. Opt. Soc. Am. B 24, 2314–2320 (2007). [CrossRef]
11. P. J. M. Vanbrabant, J. Beeckman, K. Neyts, R. James, and F. A. Fernandez, “A finite element beam propagation method for simulation of liquid crystal devices,” Opt. Express 17, 10895–10909 (2009). [CrossRef] [PubMed]
13. G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finite-element beam propagation method,” Opt. Quantum Electron. 40, 733–748 (2008). [CrossRef]
14. J. A. Fleck Jr., J. R. Morris, and M. D. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. 10, 129–160 (1976). [CrossRef]
15. B. Y. Zel’dovich and N. V. Tabiryan, “Orientational optical nonlinearity of liquid crystals,” Sov. Phys. Usp. 28, 1059–1083 (1985). [CrossRef]
16. U. A. Laudyn, K. Jaworowicz, and M. A. Karpierz, “Spatial solitons in chiral nematics,” Mol. Cryst. Liq. Cryst. 489, 214–221 (2008). [CrossRef]
17. M. Peccianti, A. Dyadyusha, M. Kaczmarek, and G. Assanto, “Tunable refraction and reflection of self-confined light beams,” Nat. Phys. 2, 737–742 (2006). [CrossRef]
19. F. A. Sala and M. A. Karpierz, “Chiral and non-chiral nematic liquid crystal reorientation induced by inhomogeneous electric fields,” J. Opt. Soc. Am. B 29 (submitted) (2012).
20. F. A. Sala and M. A. Karpierz, “Modeling of nonlinear beam propagation in chiral nematic liquid crystals,” Mol. Cryst. Liq. Cryst. 558, 176–183 (2012) [CrossRef]
21. I.-C. Khoo, Liquid Crystals (John Wiley and Sons, 2007). [CrossRef]
22. F. C. Frank, “I. Liquid crystals. On the theory of liquid crystals,” Discuss. Faraday Soc. 25, 19–28 (1958). [CrossRef]
23. C. W. Oseen, “The theory of liquid crystals,” Trans. Faraday Soc. 29, 883–899 (1933). [CrossRef]
24. W. Baran, Z. Raszewski, R. Dabrowski, and J. Kedzierski, “Some physical properties of mesogenic 4-(trans-4’-n-alkylcyclohexyl) isothiocyanatobenzenes,” Mol. Cryst. Liq. Cryst. 123, 237–245 (1985). [CrossRef]
25. R. Dabrowski, J. Dziaduszek, and T. Szczucinski, “Mesomorphic characteristics of some new homologous series with the isothiocyanato terminal group,” Mol. Cryst. Liq. Cryst. 124, 241–257 (1985). [CrossRef]
26. D. M. Young, “Iterative methods for solving partial difference equations of elliptical type,” PhD thesis, Harvard University (1950).
27. L. Hageman and D. Young, Applied Iterative Methods (Academic Press, 1981).
28. P. de Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, 1993).