**Part 2**

**Technical Schemes and Processes** 

60 Features of Liquid Crystal Display Materials and Processes

Le, H. P. (1998). Progress and trends in ink-jet printing technology. *J. Imaging Sci. Technol*.,

Lee, H.-H.; Chou, K.-S. & Huang, K.-C. (2005). Inkjet printing of nanosized silver colloids.

Lenz, P.; Fenzl, W. & Lipowsky, R. (2011). Wetting of ring-shaped surface domains.

Lim, J. A.; Kim, J.-H.; Qiu, L.; Lee W. H.; Lee, H. S.; Kwak, D. & Cho, K. (2010). Inkjet-printed

Liou, T.-M.; Chan, C.-Y.; Fu, C.-C. & Shih, K.-C. (2008). Effects of impact inertia and surface

MacFarlane, D. L.; Narayan, V.; Tatum, J. A.; Cox, W. R.; Chen, T. & Hayes, D. J. (1994). Microjet fabrication of microlens arrays. IEEE Photonics Tech. Lett., 6, 1112-1114. Mentley, D. E. (2002). State of flat-panel display technology and future trends. *Proc. of the* 

Nallani, A. K.; Chen, T.; Hayes, D. J.; Che, W.-S. & Lee, J.-B. (2006). A method for improved

Parisse, F. & Allain, C. (1997). Drying of colloidal suspension droplets: experimental study

Perelaer, J. & Schubert, U. S. (2010). *Inkjet printing and alternative sintering of narrow conductive* 

Scandurra, A.; Indelli, G. F.; Sparta, N. G.; Galliano, F.; Ravesi, S. & Pignatro S. (2010). Low-

Shaw, J. M. & Seidler, P. F. (2001). Organic electronics: introduction. *IBM J. Res. & Dev*., 45, 3-9. Sirringhaus, H.; Kawase, T.; Friend, R. H.; Shimoda, T.; Inbasekaran, M.; Wu. W. & Woo, E.

Soltman, D. & Subramanian, V. (2008). Inkjet-printed line morphologies and temperature

Szczech B. J.; Megaridis, C. M.; Gamota, D. R. & Zhang, J. (2002). Fine-line conductor

Tseng, F.-G.; Kim, C.-J.; Ho, C.-M. (2002). A high-resolution high-frequency monolithic top-

van der Vaart, N. C.; Lifka, H.; Budzelaar, F. P. M.; Rubingh, J. E. J. M.; Hoppenbrouwers, J. J.

Xu, T.; Jin, J.; Gregory, C.; Hickman, J. J. & Boland, T. (2005). Inkjet printing of viable

matrix printed polymer OLED television. *J. Soc. Inf. Display*, 13/1, 9-16. Weon, B. M. & Je, J. H. (2010). Capillary force repels coffee-ring effect. *Phys. Rev. E*, 82,

insulating polymers. *Adv. Func. Mater.*, 20, 3292-3297.

and profile renormalization. *Langmuir*, 13, 3598-3602.

Satoi, T. (2001). Color filter manufacturing apparatus. US Patent 6,331,384 B1.

plastic electronics. Surf. Interface Annal., 42, 1163-1167.

control of the coffee ring effect. *Langmuir*, 24, 2224-2231.

*Electron. Packag. Manufact.*, 25, 26-33.

*Microelectromech. Syst.* , 11, 427-436.

mammalian cells. Biomaterials, 26, 93-99.

single-droplet organic transistors based on semiconductor nanowires embedded in

characteristics on deposited polymer droplets in microcavities. *J. Microelectromech.* 

VCSEL packaging using MEMS and ink-jet technologies. *J. Lightw. Technol.*, 24,

*tracks on flexible substrates for plastic electronic applications*, INTECH, ISBN 978-953-

temperature sintered conductive silver patterns obtained by inkjet printing for

P. (2000). High-resolution inkjet printing of all-polymer transistor circuits. *Science*,

manufacturing using drop-on-demand PZT printing technology. *IEEE Trans.* 

shooting microinjector free of satellite drops – part I: concept, design, and model. *J.* 

L.; Dijksman, J. F.; Verbeek, R. G. F. A.; van Woudenberg, R.; Vossen, F. J.; Hiddink, M. G. H.; Rosink, J. J. W. M.; Bernards, T. N. M.; Giraldo, A.; Young, N. D.; Fish, D. A.; Childs, M. J.; Steer, W. A. & George, D. S. (2005). Towards large-area full-color active-

42, 49-62.

*Nanotechnology*, 16, 2436-2441.

*Europhys. Lett*., 53, 618-624.

*Syst.*, 17, 278-287.

*IEEE*, 90, 453-459.

7619-72-5, Book Chapter 16.

1504-1512.

290, 2123-2126.

015305(R) (4 pages).

**0**

**4**

*Taiwan*

**Electromagnetic Formalisms for Optical**

**Liquid-Crystal Microstructures**

I-Lin Ho and Yia-Chung Chang

**Propagation in Three-Dimensional Periodic**

*Research Center for Applied Sciences, Academia Sinica, Taipei, Taiwan 115, R.O.C.*

Nanoscale structures have achieved novel functions in liquid crystal devices such as liquid crystal displays, optical filters, optical modulators, phase conjugated systems, optical attenuators, beam amplifiers, tunable lasers, holographic data storage and even as parts for optical logic systems over the last decades (Blinov et al. (2006; 2007); Sutkowski et al. (2006)). Many theoretical works also have been reported on liquid crystal (LC) optics. Jones method (Jones (1941)) is first proposed for an easy calculation, which stratifies the media along the cell normal while remains the transverse LC orientation uniform, and hence supplies a straightforward way to analyze the forward propagation at normal incidence. This was later followed by the extended Jones method (Lien (1997)), which allows to trace the forward waves at an oblique incidence. The Berreman method (Berreman (1972)) then provides an alternative

A further step in LC optics is to consider rigorously the LC variation both along the cell normal and along a single transverse direction, leading to a two-dimensional treatment of light propagation. This step is fulfilled by implementing the finite-difference time-domain method (Kriezis et al. (2000a); Witzigmann et al. (1998)), the vector beam propagation method (Kriezisa & Elston (1999 ); Kriezis & Elston (2000b)), coupled-wave theory (Galatola et al. (1994); Rokushima & Yamakita (1983)), and an extension of the Berreman approach (Zhang & Sheng (2003)), and has proven to be successful in demonstrating the strong scattering and diffractive effects on the structures with transverse LC variation lasting over

For three-dimensional LC medium with arbitrary normal and transverse LC variations, Kriezis et al. (2002) proposed a composite scheme based on the finite-difference time-domain method and the plane-wave expansion method to evaluate the light propagation in periodic liquid-crystal microstructures. Olivero & Oldano (2003) applied numerical calculations by a standard spectral method and the finite-difference frequency-domain method for electromagnetic propagation in LC cells. Glytsis & Gaylord (1987) gave three-dimensional coupled-wave diffraction algorithms via the field decomposition into ordinary and extraordinary waves, although the transverse variation of the ordinary/extraordinary axis raises the complexity. Alternatively, this work neglects the multiple reflections and gives a coupling-matrix algorithm that is much easier to manipulate algebraically for three-dimensional LC media, yet accounts for the effects of the Fresnel refraction and

**1. Introduction**

process to include forward and backward waves.

the optical-wavelength scale.

## **Electromagnetic Formalisms for Optical Propagation in Three-Dimensional Periodic Liquid-Crystal Microstructures**

I-Lin Ho and Yia-Chung Chang

*Research Center for Applied Sciences, Academia Sinica, Taipei, Taiwan 115, R.O.C. Taiwan*

#### **1. Introduction**

Nanoscale structures have achieved novel functions in liquid crystal devices such as liquid crystal displays, optical filters, optical modulators, phase conjugated systems, optical attenuators, beam amplifiers, tunable lasers, holographic data storage and even as parts for optical logic systems over the last decades (Blinov et al. (2006; 2007); Sutkowski et al. (2006)). Many theoretical works also have been reported on liquid crystal (LC) optics. Jones method (Jones (1941)) is first proposed for an easy calculation, which stratifies the media along the cell normal while remains the transverse LC orientation uniform, and hence supplies a straightforward way to analyze the forward propagation at normal incidence. This was later followed by the extended Jones method (Lien (1997)), which allows to trace the forward waves at an oblique incidence. The Berreman method (Berreman (1972)) then provides an alternative process to include forward and backward waves.

A further step in LC optics is to consider rigorously the LC variation both along the cell normal and along a single transverse direction, leading to a two-dimensional treatment of light propagation. This step is fulfilled by implementing the finite-difference time-domain method (Kriezis et al. (2000a); Witzigmann et al. (1998)), the vector beam propagation method (Kriezisa & Elston (1999 ); Kriezis & Elston (2000b)), coupled-wave theory (Galatola et al. (1994); Rokushima & Yamakita (1983)), and an extension of the Berreman approach (Zhang & Sheng (2003)), and has proven to be successful in demonstrating the strong scattering and diffractive effects on the structures with transverse LC variation lasting over the optical-wavelength scale.

For three-dimensional LC medium with arbitrary normal and transverse LC variations, Kriezis et al. (2002) proposed a composite scheme based on the finite-difference time-domain method and the plane-wave expansion method to evaluate the light propagation in periodic liquid-crystal microstructures. Olivero & Oldano (2003) applied numerical calculations by a standard spectral method and the finite-difference frequency-domain method for electromagnetic propagation in LC cells. Glytsis & Gaylord (1987) gave three-dimensional coupled-wave diffraction algorithms via the field decomposition into ordinary and extraordinary waves, although the transverse variation of the ordinary/extraordinary axis raises the complexity. Alternatively, this work neglects the multiple reflections and gives a coupling-matrix algorithm that is much easier to manipulate algebraically for three-dimensional LC media, yet accounts for the effects of the Fresnel refraction and

in the *xz* plane at angle angle *θ* related to *z* axis, it specifies*k* = (*k*0*sinθ*, 0, *k*0*cosθ*), extended Jones Matrix can relate the electric fields at the bottom of the *th* layer to the fields at the top

<sup>65</sup> Electromagnetic Formalisms for Optical Propagation

,0

�

; **A** =

*<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o n*2 *e*

; **<sup>J</sup>** = **<sup>A</sup><sup>Ξ</sup><sup>A</sup>**−<sup>1</sup>

� *ex*<sup>1</sup> *ex*<sup>2</sup> *ey*<sup>1</sup> *ey*<sup>2</sup> �

; (3)

cos<sup>2</sup> *θ<sup>o</sup>* sin2 *φ<sup>o</sup>*

<sup>+</sup> *<sup>ε</sup>zx*� � *kx*

transforms the electric fields at the bottom of the *th*

*Ex Ey* � 0

⎤

⎤

<sup>0</sup> 2 cos *<sup>θ</sup>* cos *<sup>θ</sup>*+*np* cos *<sup>θ</sup><sup>p</sup>*

0 <sup>2</sup>*np* cos *<sup>θ</sup><sup>p</sup>*

cos *θ*+*np* cos *θ<sup>p</sup>*

*k*0 *kz*<sup>2</sup> *k*0

� *k*<sup>2</sup> *x k*2 0

− *εyzεzy* (5)

*εyz* (6)

*εzy* (7)

<sup>+</sup> *<sup>ε</sup>xz*�

⎦ (10)

⎦ (11)

�1/2

(1)

(2)

(4)

(8)

(9)

of the *th* layer of each strip by:

**Ξ** = �

*kz*<sup>1</sup> *k*0 = � *n*2 <sup>0</sup> <sup>−</sup> *<sup>k</sup>*<sup>2</sup> *x k*2 0

*kz*<sup>2</sup> *k*0

*ex*<sup>1</sup> =

*ey*<sup>1</sup> =

*ex*<sup>2</sup> =

*ey*<sup>2</sup> =

<sup>=</sup> <sup>−</sup> *<sup>ε</sup>xz εzz kx k*0 + *none εzz* �

> � *k*2 *z*1 *k*2 0 + *k*2 *x k*2 0

� *k*2 *x k*2 0

� − *k*2 *x k*2 0

� − *k*2 *z*2 *k*2 0

(1) can be understood as follow. **A**−<sup>1</sup>

� *Ex Ey* �

<sup>−</sup> *<sup>ε</sup>zz*�

<sup>+</sup> *<sup>ε</sup>zz*�

with

� *Ex Ey* �

in Three-Dimensional Periodic Liquid-Crystal Microstructures

,*dz*

0 exp (*ikz*2*dz*)

<sup>−</sup> *<sup>ε</sup>yy*� �*k*<sup>2</sup>

*εyx* +

<sup>+</sup> *<sup>ε</sup>xx*� � *<sup>k</sup>*<sup>2</sup>

*εxy* −

*x k*2 0

electric fields ( = 0) and the emitted electric fields ( = *N* + 1) is given by

⎡ ⎣

⎡ ⎣

*N*+1

**J***ent* =

**J***ext* =

*εzz* − � <sup>1</sup> <sup>−</sup> *<sup>n</sup>*<sup>2</sup>

�*kx k*0 *kz*<sup>1</sup> *k*0

> � *kx k*0 *kz*<sup>2</sup> *k*0

<sup>−</sup> *<sup>ε</sup>zz*�

*x k*2 0 <sup>−</sup> *<sup>ε</sup>zz*�

<sup>+</sup> *<sup>ε</sup>zx*�

+ � *kx k*0 *kz*<sup>2</sup> *k*0

Here, *k*<sup>0</sup> = *ω*/*c* = 2*π*/*λ* with *λ* the wavelength of the incident light in free space. *dz* is the thickness of the the *th* layer. *θ<sup>o</sup>* and *φ<sup>o</sup>* are the orientation angles of the LC director defined in the spherical coordinate. *<sup>ε</sup>i*,*j*∈{*x*,*y*,*z*} is the dielectric tensors defined in appendix A. Equation

layer into the eigen-mode fields. **Ξ** then propagates the eigen-mode fields from the bottom of the *th* layer to the top of the *th* layer through the distance *dz*. Finally, **A** transform the eigen-mode fields at the top of the the *th* layer back into the electric fields at the top of the *th* layer, which is equal to the electric fields at the bottom of the ( + 1)*th* layer by boundary condition. Grouping all layers, the extended Jones matrix formula that relates the incident

<sup>=</sup> **<sup>J</sup>***ext***J***N***J***N*−1...**J**2**J**1**J***ent* �

2 cos *θ<sup>p</sup>* cos *<sup>θ</sup>p*+*np* cos *<sup>θ</sup>* 0

2*np* cos *θ* cos *<sup>θ</sup>p*+*np* cos *<sup>θ</sup>* 0

<sup>+</sup> *<sup>ε</sup>xz*�

exp (*ikz*1*dz*) 0

�1/2

= **J** � *Ex Ey* �

the single reflection at the surfaces of the media. The detailed derivations are described in appendix A. Furthermore, analogous with the Berreman approach (Berreman (1972)) to consider the multiple reflections for one-dimensional layered media (i.e. stratifying the media along the cell normal while remaining the transverse LC orientation uniform), another supplementary formulae including the influences of multiple reflections for three-dimensional media (i.e. stratifying the media along the cell normal and simultaneously including the varying LC orientation along the transverse) are also addressed in the appendix A. The program code of wolfram mathematica for coupling-matrix method is appended in appendix B for references.

#### **2. Extended Jones matrix method revisited**

Fig. 1. (a) Schematic depiction of one unit cell of the periodic LC structures. (b) Stratification of the cell along the cell normal *z*ˆ with remaining the real transverse *x*ˆ(*y*ˆ) profile as in coupling-matrix method.(c) Decompose the cell along the transverse direction *x*ˆ(*y*ˆ) into independent strips, and treat the stratification of each stripe with uniform transverse profiles, as in (extended) Jones matrix method.

In this section, extended Jones matrix method is revisited first due to its similar underlying concepts can supply an accessibility to understand the coupling-matrix method. In the extended Jones matrix method, the liquid crystal cell (Figure 1(a)) is decomposed into multiple one-dimensional (*z*) independent stripes (Figure 1(c)), treating the transverse LC orientation uniform within each stripes and being irrelevant each other. Each stripe is further divided into *N* layers along the *z* direction, including two separate polarizer and analyzer layers. In the layer, there are four eigen-mode waves: two transmitted and two reflected waves; while at the interface of the layer, the boundary condition is that the tangential components of the electric field are continuous. Without loss of generality, considering the propagation of waves in the *xz* plane at angle angle *θ* related to *z* axis, it specifies*k* = (*k*0*sinθ*, 0, *k*0*cosθ*), extended Jones Matrix can relate the electric fields at the bottom of the *th* layer to the fields at the top of the *th* layer of each strip by:

$$\begin{bmatrix} E\_{\mathcal{X}} \\ E\_{\mathcal{Y}} \end{bmatrix}\_{\ell, dz\_{\ell}} = \mathbf{J}\_{\ell} \begin{bmatrix} E\_{\mathcal{X}} \\ E\_{\mathcal{Y}} \end{bmatrix}\_{\ell, 0}; \ \mathbf{J}\_{\ell} = \mathbf{A}\_{\ell} \boldsymbol{\Xi}\_{\ell} \mathbf{A}\_{\ell}^{-1} \tag{1}$$

with

2 Will-be-set-by-IN-TECH

the single reflection at the surfaces of the media. The detailed derivations are described in appendix A. Furthermore, analogous with the Berreman approach (Berreman (1972)) to consider the multiple reflections for one-dimensional layered media (i.e. stratifying the media along the cell normal while remaining the transverse LC orientation uniform), another supplementary formulae including the influences of multiple reflections for three-dimensional media (i.e. stratifying the media along the cell normal and simultaneously including the varying LC orientation along the transverse) are also addressed in the appendix A. The program code of wolfram mathematica for coupling-matrix method is appended in

(a) (b)

Fig. 1. (a) Schematic depiction of one unit cell of the periodic LC structures. (b) Stratification of the cell along the cell normal *z*ˆ with remaining the real transverse *x*ˆ(*y*ˆ) profile as in coupling-matrix method.(c) Decompose the cell along the transverse direction *x*ˆ(*y*ˆ) into independent strips, and treat the stratification of each stripe with uniform transverse profiles,

In this section, extended Jones matrix method is revisited first due to its similar underlying concepts can supply an accessibility to understand the coupling-matrix method. In the extended Jones matrix method, the liquid crystal cell (Figure 1(a)) is decomposed into multiple one-dimensional (*z*) independent stripes (Figure 1(c)), treating the transverse LC orientation uniform within each stripes and being irrelevant each other. Each stripe is further divided into *N* layers along the *z* direction, including two separate polarizer and analyzer layers. In the layer, there are four eigen-mode waves: two transmitted and two reflected waves; while at the interface of the layer, the boundary condition is that the tangential components of the electric field are continuous. Without loss of generality, considering the propagation of waves

Sent S1 S2

> Jent J1 J2

J ent J 1 J 2

J ext J N

Jext JN

Sext SN

(c) <sup>x</sup><sup>ˆ</sup>

appendix B for references.

**2. Extended Jones matrix method revisited**

yˆ

zˆ

as in (extended) Jones matrix method.

$$\boldsymbol{\Xi}\_{\ell} = \begin{bmatrix} \exp\left(ik\_{z1}dz\_{\ell}\right) & 0\\ 0 & \exp\left(ik\_{z2}dz\_{\ell}\right) \end{bmatrix}; \quad \mathbf{A}\_{\ell} = \begin{bmatrix} \boldsymbol{e}\_{\mathbf{x}1} \ \boldsymbol{e}\_{\mathbf{x}2} \\ \boldsymbol{e}\_{\mathbf{y}1} \ \boldsymbol{e}\_{\mathbf{y}2} \end{bmatrix} \tag{2}$$

$$\frac{k\_{z1}}{k\_0} = \left(n\_0^2 - \frac{k\_x^2}{k\_0^2}\right)^{1/2};\tag{3}$$

$$\frac{k\_{z2}}{k\_0} = -\frac{\varepsilon\_{x2}}{\varepsilon\_{zz}} \frac{k\_x}{k\_0} + \frac{n\_o n\_\varepsilon}{\varepsilon\_{zz}} \left(\varepsilon\_{zz} - \left(1 - \frac{n\_\varepsilon^2 - n\_o^2}{n\_\varepsilon^2} \cos^2 \theta\_0 \sin^2 \phi\_0 \right) \frac{k\_x^2}{k\_0^2} \right)^{1/2} \tag{4}$$

$$\varepsilon\_{\mathbf{x}1} = \left(\frac{k\_{z1}^2}{k\_0^2} + \frac{k\_{\mathbf{x}}^2}{k\_0^2} - \varepsilon\_{yy}\right) \left(\frac{k\_{\mathbf{x}}^2}{k\_0^2} - \varepsilon\_{zz}\right) - \varepsilon\_{yz}\varepsilon\_{zy} \tag{5}$$

$$e\_{y1} = \left(\frac{k\_x^2}{k\_0^2} - \varepsilon\_{zz}\right)\varepsilon\_{yx} + \left(\frac{k\_x}{k\_0}\frac{k\_{z1}}{k\_0} + \varepsilon\_{zx}\right)\varepsilon\_{yz} \tag{6}$$

$$
\varepsilon\_{\mathbf{x}2} = \left(-\frac{k\_{\mathbf{x}}^2}{k\_0^2} + \varepsilon\_{\mathbf{z}2}\right)\varepsilon\_{\mathbf{xy}} - \left(\frac{k\_{\mathbf{x}}}{k\_0}\frac{k\_{\mathbf{z}2}}{k\_0} + \varepsilon\_{\mathbf{x}2}\right)\varepsilon\_{\mathbf{zy}}\tag{7}
$$

$$\varepsilon\_{y2} = \left(-\frac{k\_{z2}^2}{k\_0^2} + \varepsilon\_{xx}\right)\left(\frac{k\_x^2}{k\_0^2} - \varepsilon\_{zz}\right) + \left(\frac{k\_x}{k\_0}\frac{k\_{z2}}{k\_0} + \varepsilon\_{xx}\right)\left(\frac{k\_x}{k\_0}\frac{k\_{z2}}{k\_0} + \varepsilon\_{xz}\right) \tag{8}$$

Here, *k*<sup>0</sup> = *ω*/*c* = 2*π*/*λ* with *λ* the wavelength of the incident light in free space. *dz* is the thickness of the the *th* layer. *θ<sup>o</sup>* and *φ<sup>o</sup>* are the orientation angles of the LC director defined in the spherical coordinate. *<sup>ε</sup>i*,*j*∈{*x*,*y*,*z*} is the dielectric tensors defined in appendix A. Equation (1) can be understood as follow. **A**−<sup>1</sup> transforms the electric fields at the bottom of the *th* layer into the eigen-mode fields. **Ξ** then propagates the eigen-mode fields from the bottom of the *th* layer to the top of the *th* layer through the distance *dz*. Finally, **A** transform the eigen-mode fields at the top of the the *th* layer back into the electric fields at the top of the *th* layer, which is equal to the electric fields at the bottom of the ( + 1)*th* layer by boundary condition. Grouping all layers, the extended Jones matrix formula that relates the incident electric fields ( = 0) and the emitted electric fields ( = *N* + 1) is given by

$$
\begin{bmatrix} E\_X \\ E\_Y \end{bmatrix}\_{N+1} = \mathbf{J}\_{ext} \mathbf{J}\_N \mathbf{J}\_{N-1} \dots \mathbf{J}\_2 \mathbf{J}\_1 \mathbf{J}\_{ent} \begin{bmatrix} E\_X \\ E\_Y \end{bmatrix}\_0 \tag{9}
$$

$$\mathbf{J}\_{ent} = \begin{bmatrix} \frac{2\cos\theta\_p}{\cos\theta\_p + n\_p\cos\theta} & 0\\ 0 & \frac{2\cos\theta}{\cos\theta + n\_p\cos\theta\_p} \end{bmatrix} \tag{10}$$

$$\mathbf{J}\_{\varepsilon xt} = \begin{bmatrix} \frac{2n\_p \cos \theta}{\cos \theta\_p + n\_p \cos \theta} & 0\\ 0 & \frac{2n\_p \cos \theta\_p}{\cos \theta + n\_p \cos \theta\_p} \end{bmatrix} \tag{11}$$

**G** =

⎡ ⎢ ⎢ ⎣

*ε*˜*xzε*˜−<sup>1</sup>

−*ε*˜*yzε*˜ −1

Jones method: the (**T**(*a*)

**T**(*a*)

as *E*<sup>+</sup>

*<sup>q</sup>* and *M* <sup>+</sup>

transform matrix **T**(*i*)

**<sup>f</sup>**ˆ*t*,0 = [*ex*,0 *hy*,0 *ey*,0 *hx*,0]

In this context, the notation

*n*˜ *xε*˜ −1

*n*˜ *<sup>y</sup>ε*˜−<sup>1</sup>

*zz ε*˜*zx n*˜ *<sup>x</sup>ε*˜−<sup>1</sup>

*zz ε*˜*zx n*˜ *<sup>y</sup>ε*˜

*zz <sup>ε</sup>*˜*zx* <sup>−</sup> *<sup>ε</sup>*˜*xx* <sup>+</sup> *<sup>n</sup>*˜ *yn*˜ *<sup>y</sup> <sup>ε</sup>*˜*xzε*˜−<sup>1</sup>

in Three-Dimensional Periodic Liquid-Crystal Microstructures

*zz ε*˜*zx* + *ε*˜*yx* + *n*˜ *xn*˜ *<sup>y</sup>* −*ε*˜*yzε*˜

indexes *<sup>α</sup>*, *<sup>β</sup>* are arranged by the relation *<sup>M</sup>* <sup>∼</sup> *<sup>ε</sup>*˜*ij*

tangential components of fields **f**ˆ*t*, = [*ex*,

the ( + 1)*th* layer as in Equation (13).

⎡ ⎢ ⎢ ⎣

*ex*,0 *hy*,0 *ey*,0 *hx*,0 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣

<sup>≡</sup> **<sup>T</sup>**(*i*) *εI*

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

Here, ˙**n***<sup>y</sup>* and ˙**n***<sup>x</sup>* are *NgNh* <sup>×</sup> *NgNh* diagonal matrices with normalized elements *nyh*

respectively. *ξ*−<sup>1</sup> is the diagonal matrix with elements 1/*ξgh* (not the inverse of the matrix *ξ*),

*<sup>q</sup>* and *M* <sup>−</sup>

*<sup>q</sup>* ( *E*− *zz <sup>n</sup>*˜ *<sup>x</sup>* <sup>−</sup> <sup>1</sup> *<sup>n</sup>*˜ *<sup>x</sup>ε*˜−<sup>1</sup>

*zz n*˜ *<sup>x</sup> n*˜ *<sup>y</sup>ε*˜

*zz <sup>n</sup>*˜ *<sup>x</sup>* <sup>−</sup>*ε*˜*yzε*˜−<sup>1</sup>

−1

*zz n*˜ *<sup>x</sup> ε*˜*xzε*˜

<sup>67</sup> Electromagnetic Formalisms for Optical Propagation

(or *Mgh*) describing the wave along **n***gh*. *n*˜ *<sup>x</sup>* (*n*˜ *<sup>y</sup>*) are *NgNh* × *NgNh* diagonal matrices with

are calculated by Equations (14-15). *<sup>ε</sup>*˜*ij*∈{*x*,*y*,*z*} are *NgNh* <sup>×</sup> *NgNh* matrices with elements *εij*,*αβ* being the Fourier transform of the spatial dielectric coefficients *εij*(*x*, *y*; *z*), in which the

derived in appendix A). Above *Ng*(*h*) define the number of considered total Fourier orders *g* (*h*) in the *x* (*y*) direction. 1 represents the *NgNh* × *NgNh* identity matrix. One may understand the Equation (17) for the *th* layer by the similar way as described in extended

*hy*, *ey*, *hx*,]

term describes eigen-mode propagation over the distance *dz* (thickness of the *th* layer); the

For the matrices **S***ent* and **S***ext* defined for the (isotropic) uniform incident ( = 0) and emitted ( = *N* + 1) regions, respectively, the eigen-modes are specially chosen (and symbolized)

the physical forward (backward) *TE* and *TM* waves as the above-mentioned. In which the

**n**˙ *y* **n**˙ *x* **n**˙ *y* **n**˙ *x* **<sup>n</sup>**˙ *<sup>y</sup>ξ εI***n**˙ *<sup>x</sup>ξ*−<sup>1</sup> <sup>−</sup>**n**˙ *<sup>y</sup><sup>ξ</sup>* <sup>−</sup>*εI***n**˙ *<sup>x</sup>ξ*−<sup>1</sup> −**n**˙ *<sup>x</sup>* **n**˙ *<sup>y</sup>* −**n**˙ *<sup>x</sup>* **n**˙ *<sup>y</sup>* **<sup>n</sup>**˙ *<sup>x</sup><sup>ξ</sup>* <sup>−</sup>*εI***n**˙ *<sup>y</sup>ξ*−<sup>1</sup> <sup>−</sup>**n**˙ *<sup>x</sup>ξ εI***n**˙ *<sup>y</sup>ξ*−<sup>1</sup>

*<sup>t</sup>* for the isotropic incident region ( = 0) is given as:

 term then is the inversely coordinate transformation from the eigen-mode components back to the spatial tangential components of fields at the next interface. Considering the continuum of tangential fields on interfaces, these fields emitted from the *th* layer hence can be straightforwardly treated as the incident fields **f**ˆ*t*,<sup>+</sup><sup>1</sup> for the ( + 1)*th* layer, and allow to follow the next transfer matrix **S**<sup>+</sup><sup>1</sup> to describe the sequential propagations of fields through

interface into the orthogonal components of the eigen-modes in the *th* layer; the *exp* �

−1

*NgNh* diagonal elements *nxg* (*nyh*) being the same (*g*, *h*) sequence as that of

−1

*zz <sup>ε</sup>*˜*zy* <sup>−</sup>*n*˜ *<sup>x</sup>ε*˜−<sup>1</sup>

*zz ε*˜*zy* − *ε*˜*xy* − *n*˜ *yn*˜ *<sup>x</sup>* −*ε*˜*xzε*˜

*zz ε*˜*zy* + *ε*˜*yy* − *n*˜ *xn*˜ *<sup>x</sup> ε*˜*yzε*˜

*zz ε*˜*zy* −*n*˜ *<sup>y</sup>ε*˜

−1

*<sup>E</sup>* (or *<sup>M</sup>* ) denotes the *NgNh* <sup>×</sup> 1 vector with components *Egh*

*E*, i.e. *Mgh* ∼ ∑*g*�

)−<sup>1</sup> term represents the coordinate transformation from the spatial

*<sup>q</sup>* ) (Ho et al. (2011); Rokushima & Yamakita (1983)), representing

⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

*<sup>ε</sup><sup>I</sup>* between the eigen-mode components and the tangential components

*zz n*˜ *<sup>y</sup>*

*E* and *M* , and

)(*h*−*h*�

⎤ ⎥ ⎥ ⎦

(18)

)*Eg*�*h*� (

*iκa dz* �

(19)

and *nxg mgh*

*mgh*

−1 *zz n*˜ *<sup>y</sup>*

−1 *zz n*˜ *<sup>y</sup>* + 1

*<sup>h</sup>*� *<sup>ε</sup>ij*,(*g*−*g*�

*<sup>t</sup>* denoted by Equations (46)-(47) at *th*

−1 *zz n*˜ *<sup>y</sup>*

with *<sup>θ</sup><sup>p</sup>* <sup>=</sup> sin−1(sin *<sup>θ</sup>*/� � *np* � ) in which � � *np* � stands for the average of the real parts of the two indices of refraction (*ne* and *no*) of the polarizer. The total transmission for the stripe is calculated by

$$trans. = \frac{|E\_{\rm x,N+1}|^2 + \cos^2 \theta \cdot |E\_{\rm y,N+1}|^2}{|E\_{\rm x,0}|^2 + \cos^2 \theta \cdot |E\_{\rm y,0}|^2} \tag{12}$$

The total transmission of the three-dimensional LC media then can be evaluated by summing up the contributions from the individual stripe.

#### **3. Coupling matrix method**

Parallel to the equation (9) by one-dimensional treatments for strips, an analogous coupling-matrix formulae for the propagations of waves through the three-dimensional periodic microstructures can be given as:

$$
\begin{bmatrix}
\vec{E}\_{q,N+1}^{+} \\
\vec{M}\_{q,N+1}^{+} \\
\vec{E}\_{q,N+1}^{-} \\
\vec{M}\_{q,N+1}^{-}
\end{bmatrix} = \mathbf{S}\_{ext}\mathbf{S}\_{N}...\mathbf{S}\_{2}\mathbf{S}\_{1}\mathbf{S}\_{ext}
\begin{bmatrix}
\vec{E}\_{q,0}^{+} \\
\vec{M}\_{q,0}^{+} \\
\vec{E}\_{q,0}^{-} \\
\vec{M}\_{q,0}^{-}
\end{bmatrix} \tag{13}
$$

Here, *E*<sup>+</sup> *<sup>q</sup>*, and *<sup>M</sup>* <sup>+</sup> *<sup>q</sup>*, ( *E*− *<sup>q</sup>*, and *<sup>M</sup>* <sup>−</sup> *q*, ) represent the physical forward (backward) *TE* and *TM* fields, i.e. transverse electric and transverse magnetic fields corresponding to the planes of the diffraction waves in the incident ( = 0) and emitted ( = *N* + 1) regions. In which the components of the vectors *E*<sup>+</sup> *<sup>q</sup>*,, *<sup>M</sup>* <sup>+</sup> *q*, , *E*− *q*, , or *M* <sup>−</sup> *<sup>q</sup>*, define the diffraction waves along the direction **<sup>n</sup>***gh* = *nxgı*ˆ+ *nyh* <sup>ˆ</sup>*<sup>j</sup>* + *<sup>ξ</sup>gh* <sup>ˆ</sup> *k*:

$$n\_{\mathcal{X}\mathcal{g}} = n\_I \sin \theta \cos \phi - g \frac{\lambda}{\Lambda\_{\mathcal{X}}} \tag{14}$$

$$n\_{yh} = n\_I \sin \theta \sin \phi - h \frac{\lambda}{\Lambda\_y} \tag{15}$$

$$\mathfrak{E}\_{\mathcal{g}h} = \sqrt{\varepsilon\_{I(E)} - n\_{\mathcal{y}h} n\_{\mathcal{y}h} - n\_{\mathcal{X}\mathcal{g}} n\_{\mathcal{X}\mathcal{g}}} \tag{16}$$

with *ε<sup>I</sup>* = *n*<sup>2</sup> *<sup>I</sup>* (*ε<sup>E</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup> *<sup>E</sup>*) being the dielectric coefficient in the incident (emitted) region. Note that the components with imaginary *ξgh* values are ignored for studied cases due to the decaying natures along the electromagnetic propagations parallel to the *z* direction. Λ*x* (Λ*y*) is the periodicity of the LC structure along the *<sup>x</sup>* (*y*) direction. **<sup>S</sup>**∈{1∼*N*} is the matrix representing the propagations of waves through the *th* structured layer. It consists of the matrix **T**(*a*) , which is the (column) eigen-vector matrix of the characteristic matrix **G** for the *th* layer , and the diagonal matrix *exp* � *iκ* (*a*) *dz* � relates to the eigen-value *κ* (*a*) of **G** with dimensionless *dz* = *dzk*0:

$$\mathbf{S}\_{\ell} = \mathbf{T}\_{\ell}^{(a)} \exp\left[i\kappa\_{\ell}^{(a)} \overline{d\boldsymbol{z}\_{\ell}}\right] (\mathbf{T}\_{\ell}^{(a)})^{-1} \tag{17}$$

4 Will-be-set-by-IN-TECH

two indices of refraction (*ne* and *no*) of the polarizer. The total transmission for the stripe is

The total transmission of the three-dimensional LC media then can be evaluated by summing

Parallel to the equation (9) by one-dimensional treatments for strips, an analogous coupling-matrix formulae for the propagations of waves through the three-dimensional

= **S***ext***S***N*...**S**2**S**1**S***ent*

fields, i.e. transverse electric and transverse magnetic fields corresponding to the planes of the diffraction waves in the incident ( = 0) and emitted ( = *N* + 1) regions. In which

*nxg* = *nI* sin *θ* cos *φ* − *g*

� *iκ* (*a*) *dz* �

 *exp* � *iκ* (*a*) *dz* � (**T**(*a*)

*nyh* <sup>=</sup> *nI* sin *<sup>θ</sup>* sin *<sup>φ</sup>* <sup>−</sup> *<sup>h</sup> <sup>λ</sup>*

Note that the components with imaginary *ξgh* values are ignored for studied cases due to the decaying natures along the electromagnetic propagations parallel to the *z* direction. Λ*x* (Λ*y*) is the periodicity of the LC structure along the *<sup>x</sup>* (*y*) direction. **<sup>S</sup>**∈{1∼*N*} is the matrix representing the propagations of waves through the *th* structured layer. It consists of the

, which is the (column) eigen-vector matrix of the characteristic matrix **G** for the

, or *M* <sup>−</sup>

*λ* Λ*x*

Λ*y*

*<sup>E</sup>*) being the dielectric coefficient in the incident (emitted) region.

<sup>2</sup> <sup>+</sup> cos<sup>2</sup> *<sup>θ</sup>* ·

<sup>2</sup> <sup>+</sup> cos<sup>2</sup> *<sup>θ</sup>* ·

� �*Ey*,*N*+<sup>1</sup> � � 2

� �*Ey*,0 � �

> ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0

) represent the physical forward (backward) *TE* and *TM*

*εI*(*E*) − *nyhnyh* − *nxgnxg* (16)

relates to the eigen-value *κ*

(*a*)

)−<sup>1</sup> (17)

of **G** with

*<sup>q</sup>*, define the diffraction waves along the

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

stands for the average of the real parts of the

<sup>2</sup> (12)

(13)

(14)

(15)

*np* �


with *<sup>θ</sup><sup>p</sup>* <sup>=</sup> sin−1(sin *<sup>θ</sup>*/� �

**3. Coupling matrix method**

*<sup>q</sup>*, and *<sup>M</sup>* <sup>+</sup>

the components of the vectors

direction **<sup>n</sup>***gh* = *nxgı*ˆ+ *nyh* <sup>ˆ</sup>*<sup>j</sup>* + *<sup>ξ</sup>gh* <sup>ˆ</sup>

*<sup>I</sup>* (*ε<sup>E</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup>

*th* layer , and the diagonal matrix *exp*

dimensionless *dz* = *dzk*0:

calculated by

Here, *E*<sup>+</sup>

with *ε<sup>I</sup>* = *n*<sup>2</sup>

matrix **T**(*a*)

*np* �

up the contributions from the individual stripe.

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,*N*+1 *M* <sup>+</sup> *q*,*N*+1 *E*− *q*,*N*+1 *M* <sup>−</sup> *q*,*N*+1

*<sup>q</sup>*, and *<sup>M</sup>* <sup>−</sup>

*q*,

*E*<sup>+</sup> *<sup>q</sup>*,, *<sup>M</sup>* <sup>+</sup> *q*, , *E*− *q*,

*k*:

*<sup>ξ</sup>gh* <sup>=</sup> �

**<sup>S</sup>** <sup>=</sup> **<sup>T</sup>**(*a*)

periodic microstructures can be given as:

*<sup>q</sup>*, ( *E*− ) in which � �

*trans*. <sup>=</sup> <sup>|</sup>*Ex*,*N*+1<sup>|</sup>

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

$$\mathbf{G}\_{\ell} = \begin{bmatrix} \boldsymbol{\tilde{n}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} & \boldsymbol{\tilde{n}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{n}}\_{\mathcal{X}} - 1 & \boldsymbol{\tilde{n}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Y}} & -\boldsymbol{\tilde{n}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{n}}\_{\mathcal{Y}} \\ \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}} - \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} + \boldsymbol{\tilde{n}}\_{\mathcal{Y}} \boldsymbol{\tilde{n}}\_{\mathcal{Y}} & \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} & \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Y}} - \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Y}} \boldsymbol{\tilde{n}}\_{\mathcal{X}} & -\boldsymbol{\tilde{\varepsilon}}\_{\mathcal{X}} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}}^{-1} \boldsymbol{\tilde{n}}\_{\mathcal{Y}} \\ \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Y}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{\varepsilon}}\_{\mathcal{Z}x} & \boldsymbol{\tilde{n}}\_{\mathcal{Y}} \boldsymbol{\tilde{\varepsilon}}\_{zz}^{-1} \boldsymbol{\tilde{n}}\_{\mathcal{X}} & \boldsymbol{\tilde{n}}\_{\mathcal{Y}} \boldsymbol{\tilde$$

In this context, the notation *<sup>E</sup>* (or *<sup>M</sup>* ) denotes the *NgNh* <sup>×</sup> 1 vector with components *Egh* (or *Mgh*) describing the wave along **n***gh*. *n*˜ *<sup>x</sup>* (*n*˜ *<sup>y</sup>*) are *NgNh* × *NgNh* diagonal matrices with *NgNh* diagonal elements *nxg* (*nyh*) being the same (*g*, *h*) sequence as that of *E* and *M* , and are calculated by Equations (14-15). *<sup>ε</sup>*˜*ij*∈{*x*,*y*,*z*} are *NgNh* <sup>×</sup> *NgNh* matrices with elements *εij*,*αβ* being the Fourier transform of the spatial dielectric coefficients *εij*(*x*, *y*; *z*), in which the indexes *<sup>α</sup>*, *<sup>β</sup>* are arranged by the relation *<sup>M</sup>* <sup>∼</sup> *<sup>ε</sup>*˜*ij <sup>E</sup>*, i.e. *Mgh* <sup>∼</sup> <sup>∑</sup>*g*�*h*� *<sup>ε</sup>ij*,(*g*−*g*� )(*h*−*h*� )*Eg*� *<sup>h</sup>*� ( derived in appendix A). Above *Ng*(*h*) define the number of considered total Fourier orders *g* (*h*) in the *x* (*y*) direction. 1 represents the *NgNh* × *NgNh* identity matrix. One may understand the Equation (17) for the *th* layer by the similar way as described in extended Jones method: the (**T**(*a*) )−<sup>1</sup> term represents the coordinate transformation from the spatial tangential components of fields **f**ˆ*t*, = [*ex*, *hy*, *ey*, *hx*,] *<sup>t</sup>* denoted by Equations (46)-(47) at *th* interface into the orthogonal components of the eigen-modes in the *th* layer; the *exp* � *iκa dz* � term describes eigen-mode propagation over the distance *dz* (thickness of the *th* layer); the **T**(*a*) term then is the inversely coordinate transformation from the eigen-mode components back to the spatial tangential components of fields at the next interface. Considering the continuum of tangential fields on interfaces, these fields emitted from the *th* layer hence can be straightforwardly treated as the incident fields **f**ˆ*t*,<sup>+</sup><sup>1</sup> for the ( + 1)*th* layer, and allow to follow the next transfer matrix **S**<sup>+</sup><sup>1</sup> to describe the sequential propagations of fields through the ( + 1)*th* layer as in Equation (13).

For the matrices **S***ent* and **S***ext* defined for the (isotropic) uniform incident ( = 0) and emitted ( = *N* + 1) regions, respectively, the eigen-modes are specially chosen (and symbolized) as *E*<sup>+</sup> *<sup>q</sup>* and *M* <sup>+</sup> *<sup>q</sup>* ( *E*− *<sup>q</sup>* and *M* <sup>−</sup> *<sup>q</sup>* ) (Ho et al. (2011); Rokushima & Yamakita (1983)), representing the physical forward (backward) *TE* and *TM* waves as the above-mentioned. In which the transform matrix **T**(*i*) *<sup>ε</sup><sup>I</sup>* between the eigen-mode components and the tangential components **<sup>f</sup>**ˆ*t*,0 = [*ex*,0 *hy*,0 *ey*,0 *hx*,0] *<sup>t</sup>* for the isotropic incident region ( = 0) is given as:

$$
\begin{aligned}
\begin{bmatrix}
\vec{\mathcal{E}}\_{x,0} \\
\dot{\vec{h}}\_{y,0} \\
\dot{\vec{e}}\_{y,0} \\
\dot{\vec{\mu}}\_{x,0}
\end{bmatrix} &= \begin{bmatrix}
\dot{\mathbf{n}}\_{y} & \dot{\mathbf{n}}\_{x} & \dot{\mathbf{n}}\_{y} & \dot{\mathbf{n}}\_{x} \\
\dot{\mathbf{n}}\_{y}\boldsymbol{\xi} & \varepsilon\_{I}\dot{\mathbf{n}}\_{x}\boldsymbol{\xi}^{x} & -\dot{\mathbf{n}}\_{y}\boldsymbol{\xi} - \varepsilon\_{I}\dot{\mathbf{n}}\_{x}\boldsymbol{\xi}^{x} \\
\dot{\mathbf{n}}\_{x}\boldsymbol{\xi} & -\varepsilon\_{I}\dot{\mathbf{n}}\_{y}\boldsymbol{\xi}^{x} & -\dot{\mathbf{n}}\_{x}\boldsymbol{\xi} & \varepsilon\_{I}\dot{\mathbf{n}}\_{y}\boldsymbol{\xi}^{x}
\end{bmatrix}^{-1} \begin{bmatrix}
\vec{\mathcal{E}}\_{q,0}^{+} \\
\dot{\mathbf{M}}\_{q,0}^{+} \\
\dot{\mathbf{E}}\_{q,0}^{-} \\
\dot{\mathbf{M}}\_{q,0}^{-}
\end{bmatrix} \\
&\equiv \mathbf{T}\_{\varepsilon\_{l}}^{(i)} \begin{bmatrix}
\vec{\mathcal{E}}\_{q,0}^{+} \\
\dot{\mathbf{M}}\_{q,0}^{+} \\
\dot{\mathbf{E}}\_{q,0}^{-} \\
\dot{\mathbf{M}}\_{q,0}^{-}
\end{bmatrix}
\tag{19}
$$

Here, ˙**n***<sup>y</sup>* and ˙**n***<sup>x</sup>* are *NgNh* <sup>×</sup> *NgNh* diagonal matrices with normalized elements *nyh mgh* and *nxg mgh* respectively. *ξ*−<sup>1</sup> is the diagonal matrix with elements 1/*ξgh* (not the inverse of the matrix *ξ*), in which *mgh* = (*nyhnyh* <sup>+</sup> *nxgnxg*)1/2, *<sup>ξ</sup>gh* = (*ε<sup>I</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and *<sup>ε</sup><sup>I</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup> *<sup>I</sup>* have been applied for the incident region. A similar transform for **f**ˆ*t*,*N*+<sup>1</sup> in the emitted region can be derived straightforwardly by replacing all the *ε<sup>I</sup>* in Equation (19) with *ε<sup>E</sup>* and can be denoted as **<sup>f</sup>**ˆ*t*,*N*+<sup>1</sup> <sup>=</sup> **<sup>T</sup>**(*i*) *<sup>ε</sup><sup>E</sup>* [ *E*<sup>+</sup> *<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>+</sup> *<sup>q</sup>*,*N*+<sup>1</sup> *E*− *<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>−</sup> *<sup>q</sup>*,*N*+1] *t* , with *<sup>ξ</sup>gh* = (*ε<sup>E</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and *ε<sup>E</sup>* = *n*<sup>2</sup> *E*.

**S***ent* is the matrix representing the light propagation from the incident region into the medium, and indicates the essential refraction and the reflection at the first interface of the medium. To consider these effects in a simple way, a virtual (isotropic) uniform layer, which has zero thickness and effective dielectric coefficient *ε<sup>a</sup>* = *n*<sup>2</sup> *avg*, e.g. *navg* = (*ne* + *no*)/2 for the liquid-crystal grating, is assumed to exist between the incident region and the 1*st* layer. **S***ent* thereby can be approximately evaluated as:

$$\mathbf{S}\_{ent} = \mathbf{T}\_{\varepsilon\_d}^{(i)} \begin{bmatrix} \mathbf{W}\_1'^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \tag{20}$$

0 1 2 3 4

in Three-Dimensional Periodic Liquid-Crystal Microstructures

0 1 2 3 4

0 order by RCWA ±1 order by RCWA 0 order by FDTD ±1 order by FDTD

<sup>1</sup> matrices of **G** in equation

010 �*<sup>t</sup>* and *M* <sup>+</sup>

⎦ (25)

<sup>⎦</sup> (26)

<sup>⎦</sup> (27)

⎦ (28)

*<sup>q</sup>*,0 =

(b) Thickness zLC

0 0.2 0.4 0.6 0.8 1

Diffraction efficiency

<sup>69</sup> Electromagnetic Formalisms for Optical Propagation

Fig. 2. Diffraction efficiency for the periodic LC structures at thickness *zLC* = 0 − 4*μm*, in which (a) the solid line indicates numerical results by RCWA with considering multiple reflections as in the appendix (equations 89-92), and are in agreement with those (dotted line) from the FDTD method, while (b) the solid line indicates numerical results by the RCWA with ignoring multiple reflections, yet accounting for the effects of the Fresnel refraction and the single reflection at the surfaces of the media as in equation (13), showing comparable

which relates to the eigen-vector matrices **T**0/**T**<sup>2</sup> for the isotropic incident/emitted layer in

(18) for the liquid-crystal film. Here, *zlc* = *zlc*/*k*<sup>0</sup> is the thickness of the liquid-crystal film. In this case, we simply choose the unit-amplitude normal TE incidence with respect

the associated ˙**n***x*, ˙**n***y*, *ε*, and *ξ* in **T**0/**T**<sup>2</sup> are referred to equations (14,15,16), and are given as:

⎡ ⎣

*<sup>x</sup>* 0 0 010 0 0 �<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*2/Λ<sup>2</sup>

Note we have used a small incident angle (*θ* = 10−5, *φ* = 0) to avoid the numerical instability at *<sup>θ</sup>* <sup>=</sup> 0. For the layer of liquid-crystal film, the associated ˜**n***x*, ˜**n***<sup>y</sup>* and *<sup>ε</sup>ij*∈{*x*,*y*,*z*} in **<sup>G</sup>** in

⎤

⎦ ; ˜**n***<sup>y</sup>* =

⎡ ⎣

000 000 000 ⎤

000 000 000

*<sup>x</sup>* 0 0 010 0 0 1/�<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*2/Λ<sup>2</sup>

*<sup>q</sup>*,0,−<sup>10</sup> *<sup>E</sup>*<sup>+</sup>

<sup>1</sup> and eigen-vector **<sup>T</sup>**(*a*)

*<sup>q</sup>*,0,00 *<sup>E</sup>*<sup>+</sup>

⎤ ⎦ ;*ε* =

*<sup>q</sup>*,0,10 �*<sup>t</sup>*

⎡ ⎣

⎤ ⎥

*x*

= �

⎤

. For the isotropic incident/emitted air layer (*ε* = 1),

100 010 001

*x*

⎤ ⎥

(*a*)

� *E*<sup>+</sup>

⎤

�<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*2/Λ<sup>2</sup>

1/�<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*2/Λ<sup>2</sup>

*λ*/Λ*x* 0 0 01 0 0 0 −*λ*/Λ*<sup>x</sup>*

⎦ ; ˙**n***<sup>y</sup>* =

*E*<sup>+</sup> *<sup>q</sup>*,0 =

000 �*<sup>t</sup>*

= �

⎡ ⎣

⎡ ⎢ ⎣

⎡ ⎢ ⎣ 0 order by RCWA ±1 order by RCWA 0 order by FDTD ±1 order by FDTD

(a) Thickness zLC

0 0.2 0.4 0.6 0.8 1

the equation (19), and the eigen-values *κ*

*<sup>q</sup>*,0,10 �*<sup>t</sup>*

**n**˙ *<sup>x</sup>* =

*ξ* =

*ξ*−<sup>1</sup> =

**n**˜ *<sup>x</sup>* =

⎡ ⎣

equation (18) are written out as below:

to the *xz* incident plane, i.e.

*<sup>q</sup>*,0,00 *<sup>M</sup>*<sup>+</sup>

*<sup>q</sup>*,0,−<sup>10</sup> *<sup>M</sup>*<sup>+</sup>

Diffraction efficiency

results.

� *M*<sup>+</sup>

$$
\begin{bmatrix}
\mathbf{W}'\_1 \ \mathbf{W}'\_2 \\
\mathbf{W}'\_3 \ \mathbf{W}'\_4
\end{bmatrix} = \begin{bmatrix}
(\mathbf{T}^{(i)}\_{\varepsilon\_d})^{-1} \mathbf{T}^{(i)}\_{\varepsilon\_I}
\end{bmatrix}^{-1} \tag{21}
$$

Here, **T**(*i*) *<sup>ε</sup> <sup>a</sup>* is formulated as equation (19) with the replacements of *ε<sup>I</sup>* by *εa*, *ξgh* = (*ε<sup>a</sup>* − *nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and *<sup>ε</sup><sup>a</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup> *avg*. Similar to the argument of **S***ent*, another virtual (isotropic) uniform layer is included between the emitted region and the *Nth* layer to consider the effects of refraction and the reflection at the last interface. Here, **S***ext* is approximated as:

$$\mathbf{S}\_{ext} = \begin{bmatrix} \mathbf{W}\_1^{\prime\prime -1} \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} (\mathbf{T}\_{\varepsilon\_d}^{(i)})^{-1} \tag{22}$$

$$
\begin{bmatrix}
\mathbf{W}\_1^{\prime\prime} \mathbf{W}\_2^{\prime\prime} \\
\mathbf{W}\_3^{\prime\prime} \mathbf{W}\_4^{\prime\prime}
\end{bmatrix} = \begin{bmatrix}
(\mathbf{T}\_{\varepsilon\_E}^{(i)})^{-1} \mathbf{T}\_{\varepsilon\_d}^{(i)}
\end{bmatrix}^{-1} \tag{23}
$$

Put everything together, and the propagation of fields through three-dimensional periodic microstructures hence can be evaluated as in Equation (13).

#### **4. Numerical analyses**

In this section, a simple case is applied to demonstrate the algorithms and is verified by finite-difference time-domain (FDTD) method. Consider a one-layer film (*N* = 1) with liquid-crystal orientation *θ<sup>o</sup>* = *πx*/Λ*<sup>x</sup>* = *λx*/2Λ*x*, *φ<sup>o</sup>* = *π*/2 . By the Fourier transform defined in equations (42-45), the non-zero Fourier components for the dielectric elements *εij*,*gh* are: *εxx*,00 = *n*<sup>2</sup> *<sup>o</sup>*, *<sup>ε</sup>yy*,00 = � *n*2 *<sup>o</sup>* + *n*<sup>2</sup> *e* � /2, *<sup>ε</sup>yy*,±<sup>10</sup> <sup>=</sup> � *n*2 *<sup>o</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *e* � /4, *<sup>ε</sup>yz*,±<sup>10</sup> = ±*<sup>i</sup>* � *n*2 *<sup>o</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *e* � /4, *εzz*,00 = � *n*2 *<sup>o</sup>* + *n*<sup>2</sup> *e* � /2, *<sup>ε</sup>zz*,±<sup>10</sup> <sup>=</sup> � *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* � /4. For simplicity, we only consider three Fourier components of fields, i.e. (*g*, *h*) = (±1, 0) and (0, 0), for this case. The corresponding transfer-matrix formula in equation (13) are then given as:

$$
\begin{bmatrix}
\vec{E}\_{q,N+1}^{+} \\
\vec{M}\_{q,N+1}^{+} \\
\vec{E}\_{q,N+1}^{-} \\
\vec{M}\_{q,N+1}^{-}
\end{bmatrix} = \mathbf{S}\_{\varepsilon\varepsilon t} \mathbf{S}\_{1} \mathbf{S}\_{\varepsilon\varepsilon t} \begin{bmatrix}
\vec{E}\_{q,0}^{+} \\
\vec{M}\_{q,0}^{+} \\
\vec{E}\_{q,0}^{-} \\
\vec{M}\_{q,0}^{-}
\end{bmatrix} \tag{24}
$$

6 Will-be-set-by-IN-TECH

applied for the incident region. A similar transform for **f**ˆ*t*,*N*+<sup>1</sup> in the emitted region can be derived straightforwardly by replacing all the *ε<sup>I</sup>* in Equation (19) with *ε<sup>E</sup>* and can be denoted

**S***ent* is the matrix representing the light propagation from the incident region into the medium, and indicates the essential refraction and the reflection at the first interface of the medium. To consider these effects in a simple way, a virtual (isotropic) uniform layer, which has zero

liquid-crystal grating, is assumed to exist between the incident region and the 1*st* layer. **S***ent*

�

*avg*. Similar to the argument of **S***ent*, another virtual (isotropic)

*<sup>ε</sup> <sup>a</sup>* )−1**T**(*i*) *εI* �−<sup>1</sup>

> � (**T**(*i*)

*n*2 *<sup>o</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *e* �

> ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

*<sup>ε</sup><sup>E</sup>* )−1**T**(*i*) *ε a* �−<sup>1</sup>

*<sup>ε</sup> <sup>a</sup>* is formulated as equation (19) with the replacements of *ε<sup>I</sup>* by *εa*, *ξgh* = (*ε<sup>a</sup>* −

uniform layer is included between the emitted region and the *Nth* layer to consider the effects

Put everything together, and the propagation of fields through three-dimensional periodic

In this section, a simple case is applied to demonstrate the algorithms and is verified by finite-difference time-domain (FDTD) method. Consider a one-layer film (*N* = 1) with liquid-crystal orientation *θ<sup>o</sup>* = *πx*/Λ*<sup>x</sup>* = *λx*/2Λ*x*, *φ<sup>o</sup>* = *π*/2 . By the Fourier transform defined in equations (42-45), the non-zero Fourier components for the dielectric elements *εij*,*gh*

/2, *<sup>ε</sup>yy*,±<sup>10</sup> <sup>=</sup> �

components of fields, i.e. (*g*, *h*) = (±1, 0) and (0, 0), for this case. The corresponding

= **S***ext***S**1**S***ent*

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ � **W**��−<sup>1</sup> <sup>1</sup> **0 0 0**

of refraction and the reflection at the last interface. Here, **S***ext* is approximated as:

**S***ext* =

**<sup>S</sup>***ent* <sup>=</sup> **<sup>T</sup>**(*i*) *ε a* � **W**�−<sup>1</sup> <sup>1</sup> **0 0 0**

*<sup>q</sup>*,*N*+1] *t* *<sup>I</sup>* have been

(20)

(21)

(23)

(24)

, with *<sup>ξ</sup>gh* = (*ε<sup>E</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and

*avg*, e.g. *navg* = (*ne* + *no*)/2 for the

*<sup>ε</sup> <sup>a</sup>* )−<sup>1</sup> (22)

/4, *<sup>ε</sup>yz*,±<sup>10</sup> = ±*<sup>i</sup>*

/4. For simplicity, we only consider three Fourier

� *n*2 *<sup>o</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *e* � /4,

in which *mgh* = (*nyhnyh* <sup>+</sup> *nxgnxg*)1/2, *<sup>ξ</sup>gh* = (*ε<sup>I</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and *<sup>ε</sup><sup>I</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup>

*<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>−</sup>

as **<sup>f</sup>**ˆ*t*,*N*+<sup>1</sup> <sup>=</sup> **<sup>T</sup>**(*i*)

*ε<sup>E</sup>* = *n*<sup>2</sup> *E*.

Here, **T**(*i*)

*<sup>ε</sup><sup>E</sup>* [ *E*<sup>+</sup>

*nyhnyh* <sup>−</sup> *nxgnxg*)1/2, and *<sup>ε</sup><sup>a</sup>* <sup>=</sup> *<sup>n</sup>*<sup>2</sup>

**4. Numerical analyses**

are: *εxx*,00 = *n*<sup>2</sup>

*n*2 *<sup>o</sup>* + *n*<sup>2</sup> *e* �

*εzz*,00 = �

*<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>+</sup>

thereby can be approximately evaluated as:

*<sup>q</sup>*,*N*+<sup>1</sup> *E*−

thickness and effective dielectric coefficient *ε<sup>a</sup>* = *n*<sup>2</sup>

� **W**� <sup>1</sup> **W**� 2

� **W**�� <sup>1</sup> **W**�� 2

microstructures hence can be evaluated as in Equation (13).

*n*2 *<sup>o</sup>* + *n*<sup>2</sup> *e* �

transfer-matrix formula in equation (13) are then given as:

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q*,*N*+1 *M* <sup>+</sup> *q*,*N*+1 *E*− *q*,*N*+1 *M* <sup>−</sup> *q*,*N*+1

*n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

/2, *<sup>ε</sup>zz*,±<sup>10</sup> <sup>=</sup> �

*<sup>o</sup>*, *<sup>ε</sup>yy*,00 = �

**W**�� <sup>3</sup> **W**�� 4 � = � (**T**(*i*)

**W**� <sup>3</sup> **W**� 4 � = � (**T**(*i*)

Fig. 2. Diffraction efficiency for the periodic LC structures at thickness *zLC* = 0 − 4*μm*, in which (a) the solid line indicates numerical results by RCWA with considering multiple reflections as in the appendix (equations 89-92), and are in agreement with those (dotted line) from the FDTD method, while (b) the solid line indicates numerical results by the RCWA with ignoring multiple reflections, yet accounting for the effects of the Fresnel refraction and the single reflection at the surfaces of the media as in equation (13), showing comparable results.

which relates to the eigen-vector matrices **T**0/**T**<sup>2</sup> for the isotropic incident/emitted layer in the equation (19), and the eigen-values *κ* (*a*) <sup>1</sup> and eigen-vector **<sup>T</sup>**(*a*) <sup>1</sup> matrices of **G** in equation (18) for the liquid-crystal film. Here, *zlc* = *zlc*/*k*<sup>0</sup> is the thickness of the liquid-crystal film. In this case, we simply choose the unit-amplitude normal TE incidence with respect to the *xz* incident plane, i.e. *E*<sup>+</sup> *<sup>q</sup>*,0 = � *E*<sup>+</sup> *<sup>q</sup>*,0,−<sup>10</sup> *<sup>E</sup>*<sup>+</sup> *<sup>q</sup>*,0,00 *<sup>E</sup>*<sup>+</sup> *<sup>q</sup>*,0,10 �*<sup>t</sup>* = � 010 �*<sup>t</sup>* and *M* <sup>+</sup> *<sup>q</sup>*,0 = � *M*<sup>+</sup> *<sup>q</sup>*,0,−<sup>10</sup> *<sup>M</sup>*<sup>+</sup> *<sup>q</sup>*,0,00 *<sup>M</sup>*<sup>+</sup> *<sup>q</sup>*,0,10 �*<sup>t</sup>* = � 000 �*<sup>t</sup>* . For the isotropic incident/emitted air layer (*ε* = 1), the associated ˙**n***x*, ˙**n***y*, *ε*, and *ξ* in **T**0/**T**<sup>2</sup> are referred to equations (14,15,16), and are given as:

$$
\dot{\mathbf{n}}\_{\mathbf{x}} = \begin{bmatrix} 1 \ 0 \ 0 \\ 0 \ 1 \ 0 \\ 0 \ 0 \ -1 \end{bmatrix}; \dot{\mathbf{n}}\_{\mathbf{y}} = \begin{bmatrix} 0 \ 0 \ 0 \\ 0 \ 0 \ 0 \\ 0 \ 0 \ 0 \end{bmatrix}; \varepsilon = \begin{bmatrix} 1 \ 0 \ 0 \\ 0 \ 1 \ 0 \\ 0 \ 0 \ 1 \end{bmatrix} \tag{25}
$$

$$
\tilde{\xi} = \begin{bmatrix}
\sqrt{1 - \lambda^2/\Lambda\_\chi^2} \, 0 & 0 \\
0 & 1 & 0 \\
0 & 0 \sqrt{1 - \lambda^2/\Lambda\_\chi^2} \\
\end{bmatrix} \tag{26}
$$

$$
\xi^{-1} = \begin{bmatrix}
1/\sqrt{1-\lambda^2/\Lambda\_\chi^2} \ 0 & 0 \\
0 & 1 & 0 \\
0 & 0 \ 1/\sqrt{1-\lambda^2/\Lambda\_\chi^2}
\end{bmatrix} \tag{27}
$$

Note we have used a small incident angle (*θ* = 10−5, *φ* = 0) to avoid the numerical instability at *<sup>θ</sup>* <sup>=</sup> 0. For the layer of liquid-crystal film, the associated ˜**n***x*, ˜**n***<sup>y</sup>* and *<sup>ε</sup>ij*∈{*x*,*y*,*z*} in **<sup>G</sup>** in equation (18) are written out as below:

$$
\tilde{\mathbf{n}}\_{\mathcal{X}} = \begin{bmatrix}
\lambda/\Lambda\_{\mathcal{X}} \ 0 & 0 \\
0 & 1 & 0 \\
0 & 0 - \lambda/\Lambda\_{\mathcal{X}}
\end{bmatrix}; \tilde{\mathbf{n}}\_{\mathcal{Y}} = \begin{bmatrix}
0 \ 0 \ 0 \\
0 \ 0 \ 0 \\
0 \ 0 \ 0
\end{bmatrix} \tag{28}
$$

Define variables *<sup>k</sup>*<sup>0</sup> = *<sup>ω</sup>*√*μ*0*ε*<sup>0</sup> = <sup>2</sup>*<sup>π</sup>*

*<sup>λ</sup>* , *<sup>Y</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup>

∇*<sup>i</sup>* = *∂*/*∂ri* = *∂*/*∂rik*<sup>0</sup> = ∇*i*/*k*0, and the equations can be derived as:

in Three-Dimensional Periodic Liquid-Crystal Microstructures

∇ × �*Y*0**<sup>E</sup>** <sup>=</sup> <sup>−</sup>*<sup>i</sup>*

∇ × �*Z*0**<sup>H</sup>** <sup>=</sup> *<sup>i</sup>ε*(*r*)

associated with the orientation of director (*θo*, *φo*):

*εxx* = *n*<sup>2</sup>

*εxy* = *εyx* =

*εxz* = *εzx* =

*εyz* = *εzy* =

between the projection of the director on the xy plane and x axis.

**A.2 Maxwell's equations in k-space descriptions**

*εyy* = *n*<sup>2</sup>

*εzz* = *n*<sup>2</sup>

*<sup>o</sup>* + � *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

*<sup>o</sup>* + � *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

*<sup>o</sup>* + � *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

with

= *i* ⎡ ⎣

> *ε* = ⎡ ⎣

� *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

� *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

� *n*2 *<sup>e</sup>* <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *o* �

where *ne* and *no* are extraordinary and ordinary indices of refraction of the birefringent liquid crystal, respectively, *θo* is the angle between the director and the z axis, and *φo* is the angle

Consider the general geometry illustrated in Figure 3 of stacked multi-layer two-dimensional periodic microstructures. To apply the rigorous coupled-wave theory to the stack, all of the layers have to define the same periodicity: Λ*x* along the x direction and Λ*y* along the y

*<sup>Z</sup>*<sup>0</sup> <sup>=</sup> � *<sup>ε</sup>*<sup>0</sup> *μ*0

<sup>71</sup> Electromagnetic Formalisms for Optical Propagation

�*Y*0**<sup>E</sup>**

Here, all the field components are assumed to have time dependence of exp (*iωt*) and are omitted everywhere. The relative permeability of the medium is assumed to be *μ* = 1. Note that *<sup>ε</sup>ij*∈{*x*,*y*,*z*} are defined as functions of position (*x*, *<sup>y</sup>*, *<sup>z</sup>*) and *<sup>ε</sup>ij* are of (*x*, *<sup>y</sup>*, *<sup>z</sup>*). *<sup>λ</sup>* is the vacuum wavelength of the incident wave. Variables *x*, *y*, *z* generally represent spatial positions while these appeared in suffix, e.g *<sup>ε</sup>ij*∈{*x*,*y*,*z*}, denote the orientations along the directions *<sup>x</sup>*ˆ, *<sup>y</sup>*ˆ, *<sup>z</sup>*ˆ. Moreover, the variable *<sup>i</sup>* is the imaginary constant number *<sup>i</sup>* <sup>=</sup> √−1 and that appeared in suffix, e.g. *dzi*, is an integer indexing number. For liquid crystals, the dielectric matrix *ε* is

> *εxx εxy εxz εyx εyy εyz εzx εzy εzz*

⎤

sin2 *θ<sup>o</sup>* cos<sup>2</sup> *φo*,

sin2 *θ<sup>o</sup>* sin2 *φo*,

sin2 *θ<sup>o</sup>* sin *φ<sup>o</sup>* cos *φo*,

sin *θo* cos *θo* cos *φo*,

sin *θo* cos *θo* sin *φo*,

*εxx* (*r*) *εxy* (*r*) *εxz* (*r*) *εyx* (*r*) *εyy* (*r*) *εyz* (*r*) *εzx* (*r*) *εzy* (*r*) *εzz* (*r*)

, *r* = *k*0*r*, *x* = *k*0*x*, *y* = *k*0*y*, *z* = *k*0*z*, and

<sup>⎦</sup> �*Y*0**<sup>E</sup>** (41)

⎦ (42)

cos<sup>2</sup> *θo*, (43)

∇ · **E** = 0 (38) ∇ · **B** = 0 (39)

⎤

�*Z*0**<sup>H</sup>** (40)

$$
\varepsilon\_{\rm{xx}} = \begin{bmatrix}
\varepsilon\_{\rm{xx,00}} \ \varepsilon\_{\rm{xx},-10} \ \varepsilon\_{\rm{xx},-20} \\
\varepsilon\_{\rm{xx,10}} \ \varepsilon\_{\rm{xx,00}} \ \varepsilon\_{\rm{xx},-10} \\
\varepsilon\_{\rm{xx,20}} \ \varepsilon\_{\rm{xx},10} \ \varepsilon\_{\rm{xx,00}} \\
\end{bmatrix} = \begin{bmatrix} n\_o^2 & 0 & 0 \\ 0 & n\_o^2 & 0 \\ 0 & 0 & n\_o^2 \end{bmatrix} \tag{29}
$$

$$
\varepsilon\_{yy} = \begin{bmatrix}
\frac{n\_o^2 + n\_e^2}{2} \cdot \frac{n\_o^2 - n\_e^2}{4} & 0 \\
\frac{n\_o^2 - n\_e^2}{4} \cdot \frac{n\_o^2 + n\_e^2}{2} \cdot \frac{n\_o^2 - n\_e^2}{4} \\
0 & \frac{n\_o^2 - n\_e^2}{4} \cdot \frac{n\_o^2 + n\_e^2}{2}
\end{bmatrix} \tag{30}
$$

$$
\varepsilon\_{zz} = \begin{bmatrix}
\frac{n\_e^2 + n\_e^2}{2} \frac{n\_e^2 - n\_o^2}{4} & 0 \\
\frac{n\_e^2 - n\_o^2}{4} \frac{n\_o^2 + n\_e^2}{2} \frac{n\_e^2 - n\_o^2}{4} \\
0 & \frac{n\_e^2 - n\_o^2}{4} \frac{n\_o^2 + n\_e^2}{2}
\end{bmatrix} \tag{31}
$$

$$
\varepsilon\_{yz} = \begin{bmatrix}
0 & \frac{-i\left(n\_o^2 - n\_\epsilon^2\right)}{4} & 0 \\
\frac{i\left(n\_o^2 - n\_\epsilon^2\right)}{4} & 0 & \frac{-i\left(n\_o^2 - n\_\epsilon^2\right)}{4} \\
0 & \frac{i\left(n\_o^2 - n\_\epsilon^2\right)}{4} & 0
\end{bmatrix} \tag{32}
$$

$$\begin{array}{ccccc} & 0 & \frac{\sqrt{1-\nu} - \varkappa\_f}{4} & 0 & \text{J} \\ \varepsilon\_{xy} = 0, & \varepsilon\_{xz} = 0 & & & \end{array} \tag{33}$$

Consequently, the eigen-values *κ* (*a*) <sup>1</sup> and eigen-vector **<sup>T</sup>**(*a*) <sup>1</sup> matrices of **G** can be numerically evaluated and a similar process for **S***ent* and **S***ext* can be followed straightforwardly. Together with all these definitions of matrixes in equation (24), the transmittance fields *E*<sup>+</sup> *<sup>q</sup>*,2 and *<sup>M</sup>* <sup>+</sup> *q*,2 then can be decided. Here, we set *λ* = 0.55*um*, Λ*<sup>x</sup>* = 2.0*um*, *no* = 1.5, and *ne* = 1.6. The numerical results of far-field diffractions for this case by RCWA ignoring the influences of the multiple reflections are shown in figure 2(b), and are in agreement with these obtained by FDTD. Besides, an alternative consideration described by the equations (89-92) in appendix A, which includes the multiple reflections, is shown in figure 2(a), and clarifies the effectiveness of the easy-manipulated algorithm in equation (13) for the three-dimensional periodic LC media.

#### **A. Derivation of the coupling matrix**

In this appendix, detailed derivations of the coupling matrix method are demonstrated for references.

#### **A.1 Maxwell's equations in spatial-space descriptions**

Without charges and currents, Maxwell's equations can be read as:

$$\nabla \cdot \mathbf{E} = 0\tag{34}$$

$$\nabla \cdot \mathbf{B} = 0\tag{35}$$

$$
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \tag{36}
$$

$$
\nabla \times \mathbf{B} = \mu \mu\_0 \varepsilon \varepsilon\_0 \frac{\partial \mathbf{E}}{\partial t} \tag{37}
$$

Define variables *<sup>k</sup>*<sup>0</sup> = *<sup>ω</sup>*√*μ*0*ε*<sup>0</sup> = <sup>2</sup>*<sup>π</sup> <sup>λ</sup>* , *<sup>Y</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> *<sup>Z</sup>*<sup>0</sup> <sup>=</sup> � *<sup>ε</sup>*<sup>0</sup> *μ*0 , *r* = *k*0*r*, *x* = *k*0*x*, *y* = *k*0*y*, *z* = *k*0*z*, and ∇*<sup>i</sup>* = *∂*/*∂ri* = *∂*/*∂rik*<sup>0</sup> = ∇*i*/*k*0, and the equations can be derived as:

$$
\overline{\nabla} \cdot \mathbf{E} = 0 \tag{38}
$$

$$
\overline{\nabla} \cdot \mathbf{B} = 0 \tag{39}
$$

$$
\boxed{\nabla \times \sqrt{\chi\_0}} \mathbf{E} = -i\sqrt{Z\_0} \mathbf{H} \tag{40}
$$

$$\begin{aligned} \overline{\nabla} \times \sqrt{Z\_0} \mathbf{H} &= i \overline{\mathbf{r}} \begin{pmatrix} \overline{\tau} \end{pmatrix} \sqrt{Y\_0} \mathbf{E} \\ &= i \begin{bmatrix} \overline{\varepsilon}\_{xx} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{xy} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{xz} \begin{pmatrix} \overline{\tau} \end{pmatrix} \\ &= i \begin{bmatrix} \overline{\varepsilon}\_{yx} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{yy} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{yz} \begin{pmatrix} \overline{\tau} \end{pmatrix} \\ \overline{\varepsilon}\_{zx} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{zy} \begin{pmatrix} \overline{\tau} \end{pmatrix} \ \overline{\varepsilon}\_{zz} \begin{pmatrix} \overline{\tau} \end{pmatrix} \end{aligned} \tag{41}$$

Here, all the field components are assumed to have time dependence of exp (*iωt*) and are omitted everywhere. The relative permeability of the medium is assumed to be *μ* = 1. Note that *<sup>ε</sup>ij*∈{*x*,*y*,*z*} are defined as functions of position (*x*, *<sup>y</sup>*, *<sup>z</sup>*) and *<sup>ε</sup>ij* are of (*x*, *<sup>y</sup>*, *<sup>z</sup>*). *<sup>λ</sup>* is the vacuum wavelength of the incident wave. Variables *x*, *y*, *z* generally represent spatial positions while these appeared in suffix, e.g *<sup>ε</sup>ij*∈{*x*,*y*,*z*}, denote the orientations along the directions *<sup>x</sup>*ˆ, *<sup>y</sup>*ˆ, *<sup>z</sup>*ˆ. Moreover, the variable *<sup>i</sup>* is the imaginary constant number *<sup>i</sup>* <sup>=</sup> √−1 and that appeared in suffix, e.g. *dzi*, is an integer indexing number. For liquid crystals, the dielectric matrix *ε* is associated with the orientation of director (*θo*, *φo*):

$$
\varepsilon = \begin{bmatrix}
\varepsilon\_{xx}\ \varepsilon\_{xy}\ \varepsilon\_{xz} \\
\varepsilon\_{yx}\ \varepsilon\_{yy}\ \varepsilon\_{yz} \\
\varepsilon\_{zx}\ \varepsilon\_{zy}\ \varepsilon\_{zz}
\end{bmatrix} \tag{42}
$$

with

8 Will-be-set-by-IN-TECH

*n*2 *o*−*n*<sup>2</sup> *e* 4

*n*2 *o*+*n*<sup>2</sup> *e* 2

*n*2 *e*−*n*<sup>2</sup> *o* 4

*n*2 *o*+*n*<sup>2</sup> *e* 2

<sup>4</sup> <sup>0</sup> <sup>−</sup>*i*(*n*<sup>2</sup>

*o*−*n*<sup>2</sup> *e*) <sup>4</sup> 0

<sup>1</sup> and eigen-vector **<sup>T</sup>**(*a*)

evaluated and a similar process for **S***ent* and **S***ext* can be followed straightforwardly. Together

then can be decided. Here, we set *λ* = 0.55*um*, Λ*<sup>x</sup>* = 2.0*um*, *no* = 1.5, and *ne* = 1.6. The numerical results of far-field diffractions for this case by RCWA ignoring the influences of the multiple reflections are shown in figure 2(b), and are in agreement with these obtained by FDTD. Besides, an alternative consideration described by the equations (89-92) in appendix A, which includes the multiple reflections, is shown in figure 2(a), and clarifies the effectiveness of the easy-manipulated algorithm in equation (13) for the three-dimensional periodic LC

In this appendix, detailed derivations of the coupling matrix method are demonstrated for

∇ × **<sup>E</sup>** <sup>=</sup> <sup>−</sup> *<sup>∂</sup>***<sup>B</sup>**

∇ × **B** = *μμ*0*εε*<sup>0</sup>

*∂***E**

*o*−*n*<sup>2</sup> *e*) <sup>4</sup> 0

⎤ ⎥ <sup>⎦</sup> <sup>=</sup>

⎤ ⎥ ⎥ ⎥ ⎦

⎤ ⎥ ⎥ ⎥ ⎦

> *o*−*n*<sup>2</sup> *e*) 4

*εxy* = 0, *εxz* = 0 (33)

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

∇ · **E** = 0 (34) ∇ · **B** = 0 (35)

*<sup>∂</sup><sup>t</sup>* (36)

*<sup>∂</sup><sup>t</sup>* (37)

<sup>1</sup> matrices of **G** can be numerically

*E*<sup>+</sup>

*<sup>q</sup>*,2 and *<sup>M</sup>* <sup>+</sup>

⎡ ⎢ ⎢ ⎣ *n*2 *<sup>o</sup>* 0 0 0 *n*<sup>2</sup> *<sup>o</sup>* 0 0 0 *n*<sup>2</sup> *o* ⎤ ⎥ ⎥ ⎦

(29)

(30)

(31)

(32)

*q*,2

*εxx*,00 *εxx*,−<sup>10</sup> *εxx*,−<sup>20</sup> *εxx*,10 *εxx*,00 *εxx*,−<sup>10</sup> *εxx*,20 *εxx*,10 *εxx*,00

> *n*2 *o*−*n*<sup>2</sup> *e* <sup>4</sup> 0

*n*2 *o*+*n*<sup>2</sup> *e* 2

*n*2 *e*−*n*<sup>2</sup> *o* <sup>4</sup> 0

*n*2 *o*+*n*<sup>2</sup> *e* 2

<sup>0</sup> <sup>−</sup>*i*(*n*<sup>2</sup>

<sup>0</sup> *<sup>i</sup>*(*n*<sup>2</sup>

with all these definitions of matrixes in equation (24), the transmittance fields

<sup>0</sup> *<sup>n</sup>*<sup>2</sup> *o*−*n*<sup>2</sup> *e* 4

<sup>0</sup> *<sup>n</sup>*<sup>2</sup> *e*−*n*<sup>2</sup> *o* 4

*εxx* =

*εyy* =

*εzz* =

*εyz* =

Consequently, the eigen-values *κ*

**A. Derivation of the coupling matrix**

**A.1 Maxwell's equations in spatial-space descriptions**

Without charges and currents, Maxwell's equations can be read as:

media.

references.

⎡ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *n*2 *o*+*n*<sup>2</sup> *e* 2

*n*2 *o*−*n*<sup>2</sup> *e* 4

*n*2 *o*+*n*<sup>2</sup> *e* 2

*n*2 *e*−*n*<sup>2</sup> *o* 4

*i*(*n*<sup>2</sup> *o*−*n*<sup>2</sup> *e*)

(*a*)

$$\begin{aligned} \varepsilon\_{xx} &= n\_o^2 + \left(n\_e^2 - n\_o^2\right) \sin^2\theta\_o \cos^2\phi\_{o\prime} \\ \varepsilon\_{xy} &= \varepsilon\_{yx} = \left(n\_e^2 - n\_o^2\right) \sin^2\theta\_o \sin\phi\_o \cos\phi\_{o\prime} \\ \varepsilon\_{xz} &= \varepsilon\_{zx} = \left(n\_e^2 - n\_o^2\right) \sin\theta\_o \cos\theta\_o \cos\phi\_{o\prime} \\ \varepsilon\_{yy} &= n\_o^2 + \left(n\_e^2 - n\_o^2\right) \sin^2\theta\_o \sin^2\phi\_{o\prime} \\ \varepsilon\_{yz} &= \varepsilon\_{zy} = \left(n\_e^2 - n\_o^2\right) \sin\theta\_o \cos\theta\_o \sin\phi\_{o\prime} \\ \varepsilon\_{zz} &= n\_o^2 + \left(n\_e^2 - n\_o^2\right) \cos^2\theta\_o \end{aligned} \tag{43}$$

where *ne* and *no* are extraordinary and ordinary indices of refraction of the birefringent liquid crystal, respectively, *θo* is the angle between the director and the z axis, and *φo* is the angle between the projection of the director on the xy plane and x axis.

#### **A.2 Maxwell's equations in k-space descriptions**

Consider the general geometry illustrated in Figure 3 of stacked multi-layer two-dimensional periodic microstructures. To apply the rigorous coupled-wave theory to the stack, all of the layers have to define the same periodicity: Λ*x* along the x direction and Λ*y* along the y

with the Ψ angle between the electric field vector and the incident plane.

**f**ˆ*<sup>t</sup>* =

�*Y*0*Ey* <sup>−</sup> *<sup>∂</sup><sup>y</sup>*

�*Z*0*Hz***z**<sup>ˆ</sup> <sup>−</sup> *<sup>i</sup>*

�*Z*0*Hy* <sup>−</sup> *<sup>∂</sup><sup>y</sup>*

�*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>z</sup>* **<sup>z</sup>**<sup>ˆ</sup> <sup>+</sup> *<sup>i</sup>*

**f**ˆ*<sup>t</sup>* and **f***n*ˆ. For the component *hz*,*gh* (*z*), the equation can be derived as:

*gh*

= −*i*

*gh*

− ∑ *gh*

= −*i*

− ∑ *gh*

�*Z*0*Hx* <sup>−</sup> *<sup>∂</sup><sup>x</sup>*

�*Y*0*Ex* <sup>−</sup> *<sup>∂</sup><sup>x</sup>*

⎡ ⎢ ⎢ ⎣

*ex hy ey hx* ⎤ ⎥ ⎥ <sup>⎦</sup> , **<sup>f</sup>***n*<sup>ˆ</sup> <sup>=</sup>

at the interfaces as

spatial *x*, *y*, *z* components:

*∂x*

It is simplified to be:

*∂y*

For the component *<sup>∂</sup>ey*,*gh*(*<sup>z</sup>* )

�*Y*0*Ez* <sup>−</sup> *<sup>∂</sup><sup>z</sup>*

�*Y*0*Ey* <sup>−</sup> *<sup>∂</sup><sup>y</sup>*

∇ × �*Y*0**<sup>E</sup>** <sup>=</sup>

∇ × �*Z*0**<sup>H</sup>** <sup>=</sup>

� *∂x*

in Three-Dimensional Periodic Liquid-Crystal Microstructures

+ � *∂z*

= −*i*

� *∂x*

+ � *∂z*

�*Y*0*Ex* <sup>=</sup> ∑

�*Y*0*Ey* <sup>=</sup> ∑

= *i*

Now we can express Maxwell's equations by the (*g*, *h*) Fourier components in k-space descriptions. For simplicity, we introduce the definitions of the tangential and normal fields

<sup>73</sup> Electromagnetic Formalisms for Optical Propagation

Here,*ei*∈{*x*,*y*,*z*} <sup>=</sup>*ei* (*<sup>z</sup>*) and*hi*∈{*x*,*y*,*z*} <sup>=</sup>*hi* (*<sup>z</sup>*) are column matrices with Fourier components *ei*,*gh* (*z*) and *hi*,*gh* (*z*), respectively. In the following context, a straightforward calculation to obtain the infinite set of coupled-wave equations corresponding to the infinite Fourier components is fulfilled. First, we express Maxwell's curl equations (40)-(41) in terms of the

> �*Y*0*Ex* � **z**ˆ + � *∂y*

> > �*Y*0*Ez* � **y**ˆ

�*Z*0*Hx***x**<sup>ˆ</sup> <sup>−</sup> *<sup>i</sup>*

�*Z*0*Hz* � **y**ˆ

�*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>x</sup>* **<sup>x</sup>**<sup>ˆ</sup> <sup>+</sup> *<sup>i</sup>*

−*i* �

> −*i* �

*hz*,*gh* (*<sup>z</sup>*) exp �

*hz*,*gh* (*z*) = *nxgey*,*gh* (*z*) − *nyhex*,*gh* (*z*) (56)

−*i* �

*hx*,*gh* (*<sup>z</sup>*) exp �

�*Z*0*Hx* � **z**ˆ + � *∂y*

Next, we introduce the Fourier representations of **E**, **H**, and *ε*(*r*) as defined in Equations (44)-(47). Maxwell's curl equations (53)-(54) can thereby be regrouped by the components of

<sup>−</sup>*inxgey*,*gh* (*<sup>z</sup>*) exp �

�*Z*0*Hz* <sup>=</sup> <sup>−</sup>*i*∑

*<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

<sup>−</sup>*inyhez*,*gh* (*<sup>z</sup>*) exp �

*∂ey*,*gh* (*z*) *∂z*

�*Z*0*Hx* <sup>=</sup> <sup>−</sup>*i*∑

<sup>−</sup>*inyhex*,*gh* (*<sup>z</sup>*) exp �

*gh*

exp � −*i* �

*gh*

� *ez hz* �

�*Y*0*Ez* <sup>−</sup> *<sup>∂</sup><sup>z</sup>*

�*Z*0*Hz* <sup>−</sup> *<sup>∂</sup><sup>z</sup>*

*nxgx* + *nyhy*

*nxgx* + *nyhy*

−*i* �

*nxgx* + *nyhy*

−*i* �

*nxgx* + *nyhy*

�*Y*0*Ey* � **x**ˆ

�*Z*0*Hy***y**<sup>ˆ</sup> (53)

��

��

��

��

*nxgx* + *nyhy*

*nxgx* + *nyhy*

�� (55)

�� (57)

�*Z*0*Hy* � **x**ˆ

�*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>y</sup>* **<sup>y</sup>**<sup>ˆ</sup> (54)

(52)

Fig. 3. Geometry of three-dimensional RCWA algorithm for a multi-layer stack with two-dimensional periodic microstructures in arbitrary isotropic and birefringent material arrangement.

direction. The thickness for the *th* layer is *dz*, and these layers contribute to a total thickness of the stack *ZN* = ∑*<sup>N</sup>* <sup>=</sup><sup>1</sup> *dz*. The periodic permittivity of an individual layer in the stack can be expanded in Fourier series of the spatial harmonics as:

$$\overline{\varepsilon}\_{\overline{i}\overline{j}}\left(\overline{\mathbf{x}}, \overline{y}; \overline{z}\_{\ell}\right) = \sum\_{\mathbf{g}, \hbar} \overline{\varepsilon}\_{\overline{i}\overline{j}, \mathbf{g}\hbar}\left(\overline{z}\_{\ell}\right) \exp\left(i\frac{\mathbf{g}\lambda\overline{\mathbf{x}}}{\Lambda\_{\mathbf{x}}} + i\frac{\hbar\lambda\overline{y}}{\Lambda\_{\mathbf{y}}}\right) \tag{44}$$

$$\overline{\pi}\_{ij,gh}\left(\overline{\pi}\_{\ell}\right) = \frac{\lambda}{2\pi\Lambda\_{\overline{\lambda}}}\frac{\lambda}{2\pi\Lambda\_{\overline{\lambda}}}\int\_{0}^{\frac{2\pi\Lambda\_{\overline{\lambda}}}{\lambda}}\int\_{0}^{\frac{2\pi\Lambda\_{\overline{\lambda}}}{\lambda}}\overline{\pi}\_{ij}\left(\overline{\pi},\overline{\eta};\overline{\pi}\_{\ell}\right)\exp\left(-i\frac{g\lambda\overline{\pi}}{\Lambda\_{\overline{\lambda}}} - i\frac{h\lambda\overline{\eta}}{\Lambda\_{\overline{\lambda}}}\right)d\overline{\pi}d\overline{\eta} \tag{45}$$

A similar transform for the fields in the stack can be expressed in terms of Rayleigh expansions:

$$\sqrt{\chi\_0} \mathbf{E}\left(\overline{\mathbf{x}}, \overline{\mathbf{y}}; \overline{\mathbf{z}}\_\ell\right) = \sum\_{\mathbf{g}, \mathbf{h}} \mathbf{e}\_{\mathbf{g} \mid \mathbf{h}}\left(\overline{\mathbf{z}}\_\ell\right) \exp\left[-i\left(n\_{\mathbf{x}\mathbf{g}}\overline{\mathbf{x}} + n\_{\mathbf{y}\mathbf{h}}\overline{\mathbf{y}}\right)\right] \tag{46}$$

$$\sqrt{\mathbb{Z}\_{0}}\mathbf{H}\left(\overline{\mathbf{x}},\overline{\mathbf{y}};\overline{\mathbf{z}}\_{\ell}\right) = \sum\_{\mathcal{S},h} \mathbf{h}\_{\mathcal{g}h}\left(\overline{\mathbf{z}}\_{\ell}\right) \exp\left[-i\left(n\_{\mathcal{X}}\overline{\mathbf{x}} + n\_{\mathcal{y}h}\overline{\mathbf{y}}\right)\right] \tag{47}$$

$$n\_{\rm xg} = n\_I \sin \theta \cos \phi - g \frac{\lambda}{\Lambda\_\chi} \tag{48}$$

$$n\_{yh} = n\_I \sin \theta \sin \phi - h \frac{\lambda}{\Lambda\_Y} \tag{49}$$

where *nI* (*nE*) is the refraction index for the isotropic incident (emitted) region. *θ*, *φ* are the incident angles defined in sphere coordinates, and *z* is the normal direction for the *xy* plane of periodic structures. Here, the electric field of an incident unit-amplitude wave has been introduced by **E***inc* = **u** × *exp* (−*i***k** · **r**) as illustrated in figure 3, in which the wave vector **k** as well as the unit polarization vector **u** are given by:

$$\mathbf{k} = k\_0 n\_I \left( \sin \theta \cos \phi \mathbf{\hat{x}} + \sin \theta \sin \phi \mathbf{\hat{y}} + \cos \theta \mathbf{\hat{z}} \right) \tag{50}$$

$$\mathbf{u} = u\_x \mathbf{\hat{x}} + u\_y \mathbf{\hat{y}} + u\_z \mathbf{\hat{z}} = \left(\cos \mathbf{\hat{y}} \cos \theta \cos \phi - \sin \mathbf{\hat{y}} \sin \phi\right) \mathbf{\hat{z}}\tag{51}$$

$$+ \left( \cos \Psi \cos \theta \sin \phi + \sin \Psi \cos \phi \right) \hat{y} - \left( \cos \Psi \sin \theta \right) \hat{z}$$

10 Will-be-set-by-IN-TECH

φ

y

Λ

direction. The thickness for the *th* layer is *dz*, and these layers contribute to a total thickness

 <sup>2</sup>*π*Λ*<sup>y</sup> λ* 0

A similar transform for the fields in the stack can be expressed in terms of Rayleigh

**e***gh* (*z*) exp

**h***gh* (*z*) exp

<sup>=</sup><sup>1</sup> *dz*. The periodic permittivity of an individual layer in the stack can

*εij* (*x*, *y*; *z*) exp

 −*i* 

 −*i* 

> *λ* Λ*x*

Λ*y*

**k** = *k*0*nI* (sin *θ* cos *φx*ˆ + sin *θ* sin *φy*ˆ + cos *θz*ˆ) (50) **u** = *uxx*ˆ + *uyy*ˆ + *uzz*ˆ = (cos Ψ cos *θ* cos *φ* − sin Ψ sin *φ*) *x*ˆ (51)

 −*i gλx* Λ*x* − *i hλy* Λ*y*

*nxgx* + *nyhy*

*nxgx* + *nyhy*

Fig. 3. Geometry of three-dimensional RCWA algorithm for a multi-layer stack with two-dimensional periodic microstructures in arbitrary isotropic and birefringent material

> *i gλx* Λ*x* + *i hλy* Λ*y*

 <sup>2</sup>*π*Λ*<sup>x</sup> <sup>λ</sup>* 0

*g*,*h*

*g*,*h*

*nxg* = *nI* sin *θ* cos *φ* − *g*

*nyh* <sup>=</sup> *nI* sin *<sup>θ</sup>* sin *<sup>φ</sup>* <sup>−</sup> *<sup>h</sup> <sup>λ</sup>*

where *nI* (*nE*) is the refraction index for the isotropic incident (emitted) region. *θ*, *φ* are the incident angles defined in sphere coordinates, and *z* is the normal direction for the *xy* plane of periodic structures. Here, the electric field of an incident unit-amplitude wave has been introduced by **E***inc* = **u** × *exp* (−*i***k** · **r**) as illustrated in figure 3, in which the wave vector **k** as

+ (cos Ψ cos *θ* sin *φ* + sin Ψ cos *φ*) *y*ˆ − (cos Ψ sin *θ*) *z*ˆ

x y

k

θ ψ

u

x

be expanded in Fourier series of the spatial harmonics as:

*εij*,*gh* (*z*) exp

*λ* 2*π*Λ*y*

*<sup>Y</sup>*0**<sup>E</sup>** (*x*, *<sup>y</sup>*; *<sup>z</sup>*) <sup>=</sup> ∑

*<sup>Z</sup>*0**<sup>H</sup>** (*x*, *<sup>y</sup>*; *<sup>z</sup>*) <sup>=</sup> ∑

Λ

z z ... z

2 3

N-1

arrangement.

expansions:

of the stack *ZN* = ∑*<sup>N</sup>*

*εij* (*x*, *y*; *z*) = ∑

*<sup>ε</sup>ij*,*gh* (*<sup>z</sup>*) <sup>=</sup> *<sup>λ</sup>*

*g*,*h*

2*π*Λ*x*

well as the unit polarization vector **u** are given by:

z

Polarizer

periodic structures with

isotropic/birefringent materials

(44)

(46)

(47)

(48)

(49)

*dxdy* (45)

Polarizer

with the Ψ angle between the electric field vector and the incident plane.

Now we can express Maxwell's equations by the (*g*, *h*) Fourier components in k-space descriptions. For simplicity, we introduce the definitions of the tangential and normal fields at the interfaces as

$$\mathbf{f}\_{\vec{l}} = \begin{bmatrix} \vec{\vec{e}}\_{\mathcal{X}} \\ \vec{\vec{l}}\_{\mathcal{Y}} \\ \vec{\vec{e}}\_{\mathcal{Y}} \\ \vec{\vec{h}}\_{\mathcal{X}} \end{bmatrix}, \quad \mathbf{f}\_{\mathcal{H}} = \begin{bmatrix} \vec{\vec{e}}\_{z} \\ \vec{\vec{h}}\_{z} \end{bmatrix} \tag{52}$$

Here,*ei*∈{*x*,*y*,*z*} <sup>=</sup>*ei* (*<sup>z</sup>*) and*hi*∈{*x*,*y*,*z*} <sup>=</sup>*hi* (*<sup>z</sup>*) are column matrices with Fourier components *ei*,*gh* (*z*) and *hi*,*gh* (*z*), respectively. In the following context, a straightforward calculation to obtain the infinite set of coupled-wave equations corresponding to the infinite Fourier components is fulfilled. First, we express Maxwell's curl equations (40)-(41) in terms of the spatial *x*, *y*, *z* components:

∇ × �*Y*0**<sup>E</sup>** <sup>=</sup> � *∂x* �*Y*0*Ey* <sup>−</sup> *<sup>∂</sup><sup>y</sup>* �*Y*0*Ex* � **z**ˆ + � *∂y* �*Y*0*Ez* <sup>−</sup> *<sup>∂</sup><sup>z</sup>* �*Y*0*Ey* � **x**ˆ + � *∂z* �*Y*0*Ex* <sup>−</sup> *<sup>∂</sup><sup>x</sup>* �*Y*0*Ez* � **y**ˆ = −*i* �*Z*0*Hz***z**<sup>ˆ</sup> <sup>−</sup> *<sup>i</sup>* �*Z*0*Hx***x**<sup>ˆ</sup> <sup>−</sup> *<sup>i</sup>* �*Z*0*Hy***y**<sup>ˆ</sup> (53) ∇ × �*Z*0**<sup>H</sup>** <sup>=</sup> � *∂x* �*Z*0*Hy* <sup>−</sup> *<sup>∂</sup><sup>y</sup>* �*Z*0*Hx* � **z**ˆ + � *∂y* �*Z*0*Hz* <sup>−</sup> *<sup>∂</sup><sup>z</sup>* �*Z*0*Hy* � **x**ˆ + � *∂z* �*Z*0*Hx* <sup>−</sup> *<sup>∂</sup><sup>x</sup>* �*Z*0*Hz* � **y**ˆ = *i* �*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>z</sup>* **<sup>z</sup>**<sup>ˆ</sup> <sup>+</sup> *<sup>i</sup>* �*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>x</sup>* **<sup>x</sup>**<sup>ˆ</sup> <sup>+</sup> *<sup>i</sup>* �*Y*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>y</sup>* **<sup>y</sup>**<sup>ˆ</sup> (54)

Next, we introduce the Fourier representations of **E**, **H**, and *ε*(*r*) as defined in Equations (44)-(47). Maxwell's curl equations (53)-(54) can thereby be regrouped by the components of **f**ˆ*<sup>t</sup>* and **f***n*ˆ. For the component *hz*,*gh* (*z*), the equation can be derived as:

$$\partial\_{\overline{\pi}}\sqrt{Y\_0}E\_{\overline{\mathcal{Y}}} - \partial\_{\overline{\mathcal{Y}}}\sqrt{Y\_0}E\_{\overline{\mathcal{X}}} = \sum\_{gh} - in\_{\overline{\mathcal{X}}}e\_{y,gh}\left(\Xi\_{\ell}\right)\exp\left[-i\left(n\_{xg}\overline{\mathfrak{x}} + n\_{gh}\overline{\mathfrak{y}}\right)\right]$$

$$-\sum\_{gh} - in\_{yh}e\_{x,gh}\left(\Xi\_{\ell}\right)\exp\left[-i\left(n\_{xg}\overline{\mathfrak{x}} + n\_{gh}\overline{\mathfrak{y}}\right)\right]$$

$$= -i\sqrt{Z\_0}H\_{\overline{\mathcal{Z}}} = -i\sum\_{gh}h\_{z,gh}\left(\Xi\_{\ell}\right)\exp\left[-i\left(n\_{xg}\overline{\mathfrak{x}} + n\_{gh}\overline{\mathfrak{y}}\right)\right] \tag{55}$$

It is simplified to be:

$$h\_{\overline{\mathsf{z}}, \mathsf{g}^{\mathsf{h}}}(\overline{\mathsf{z}}\_{\ell}) = n\_{\mathsf{x}\mathcal{G}} e\_{\mathsf{y}, \mathsf{g}^{\mathsf{h}}}(\overline{\mathsf{z}}\_{\ell}) - n\_{\mathsf{y}\mathcal{H}} e\_{\mathsf{x}, \mathsf{g}^{\mathsf{h}}}(\overline{\mathsf{z}}\_{\ell}) \tag{56}$$

For the component *<sup>∂</sup>ey*,*gh*(*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

$$\partial\_{\overline{\mathcal{Y}}} \sqrt{Y\_0} E\_z - \partial\_{\overline{\mathcal{Z}}} \sqrt{Y\_0} E\_{\overline{\mathcal{Y}}} = \sum\_{g^{\text{th}}} -i n\_{g\text{fl}} e\_{z,g\text{fl}} \left( \overline{z}\_{\ell} \right) \exp \left[ -i \left( n\_{x\overline{\mathcal{X}}} \overline{x} + n\_{y\text{fl}} \overline{y} \right) \right]$$

$$- \sum\_{g\text{h}} \frac{\partial e\_{y,g\text{h}} \left( \overline{z}\_{\ell} \right)}{\partial \overline{z}} \exp \left[ -i \left( n\_{x\overline{\mathcal{X}}} \overline{x} + n\_{y\text{fl}} \overline{y} \right) \right]$$

$$= -i \sqrt{Z\_0} H\_{\overline{\mathcal{Y}}} = -i \sum\_{g\text{h}} h\_{x,g\text{h}} \left( \overline{z}\_{\ell} \right) \exp \left[ -i \left( n\_{x\overline{\mathcal{Y}}} \overline{x} + n\_{y\text{fl}} \overline{y} \right) \right] \tag{57}$$

For the component *<sup>∂</sup>hy*,*gh*(*<sup>z</sup>* )

and is simplified to be:

*∂hy*,*gh* (*z*) *∂z*

For the component *<sup>∂</sup>hx*,*gh* (*<sup>z</sup>* )

*<sup>Z</sup>*0*Hx* <sup>−</sup> *<sup>∂</sup><sup>x</sup>*

*∂z*

and is simplified to be:

*∂hx*,*gh* (*z*) *∂z*

*<sup>Z</sup>*0*Hz* <sup>−</sup> *<sup>∂</sup><sup>z</sup>*

*<sup>Z</sup>*0*Hy* <sup>=</sup> ∑

in Three-Dimensional Periodic Liquid-Crystal Microstructures

*gh*

− ∑ *m*

= *i*

= *i* ∑ *ghuv*

> +*i* ∑ *ghuv*

> +*i* ∑ *ghuv*

= −*inyhhz*,*gh* (*<sup>z</sup>*) − *<sup>i</sup>* ∑

*<sup>ε</sup>xy*,(*g*−*u*�

*gh*

= *i*

= *i* ∑ *ghuv*

> +*i* ∑ *ghuv*

> +*i* ∑ *ghuv*

= −*inxhhz*,*gh* (*<sup>z</sup>*) + *<sup>i</sup>* ∑

*<sup>ε</sup>yy*,(*g*−*u*�

+*i* ∑ *u*�*v*� − ∑ *gh*

−*<sup>i</sup>* ∑ *u*� *v*�

*<sup>Z</sup>*0*Hz* <sup>=</sup> ∑

*∂y*

*<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

<sup>−</sup>*inyhhz*,*gh* (*<sup>z</sup>*) exp

<sup>75</sup> Electromagnetic Formalisms for Optical Propagation

*<sup>ε</sup>xx*,*uvex*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>xy*,*uvey*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>xz*,*uvez*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>xx*,(*g*−*u*�

exp −*i* 

<sup>−</sup>*inxhhz*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>yx*,*uvex*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>yy*,*uvey*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>yz*,*uvez*,*gh* (*<sup>z</sup>*) exp

*<sup>ε</sup>yx*,(*g*−*u*�

)(*h*−*v*�

*<sup>v</sup>*� (*z*) + *i* ∑

*u*�*v*�

)*ey*,*u*�

)(*h*−*v*�

)(*h*−*v*�

*u*� *v*�

*<sup>v</sup>*� (*<sup>z</sup>*) − *<sup>i</sup>* ∑

*∂hy*,*gh* (*z*) *∂z*

*<sup>Y</sup>*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>x</sup>*

*u*� *v*�

)*ey*,*u*�

)(*h*−*v*�

*<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

*∂hx*,*gh* (*z*) *∂z*

*<sup>Y</sup>*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>y</sup>*

−*i* 

> −*i*

> > −*i*

−*i* 

)*ex*,*u*�

*<sup>v</sup>*� (*z*)

)(*h*−*v*�

*nxgx* + *nyhy*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

)*ez*,*u*�

*<sup>ε</sup>xz*,(*g*−*u*�

*nxgx* + *nyhy*

−*i* 

−*i* 

> −*i*

−*i* 

)*ex*,*u*�

*u*�*v*�

*<sup>v</sup>*� (*z*)

)(*h*−*v*�

)*ez*,*u*�

*<sup>ε</sup>yz*,(*g*−*u*�

exp −*i*  *nxgx* + *nyhy*

*nxgx* + *nyhy*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

(63)

*<sup>v</sup>*� (*z*) (64)

*<sup>v</sup>*� (*z*) (66)

(65)

and is simplified to be:

$$\frac{\partial e\_{y,\mathcal{S}h}(\Xi\_{\ell})}{\partial \Xi} = -i\eta\_{yh}e\_{z,\mathcal{S}h}(\Xi\_{\ell}) + ih\_{x,\mathcal{S}h}(\Xi\_{\ell}) \tag{58}$$

For the component *<sup>∂</sup>ex*,*gh* (*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

$$\begin{split} \partial\_{\overline{\pi}} \sqrt{Y\_0} E\_{\overline{x}} - \partial\_{\overline{\pi}} \sqrt{Y\_0} E\_{\overline{z}} &= \sum\_{m} \frac{\partial e\_{x, \mathcal{S}h} \left( \overline{z}\_{\ell} \right)}{\partial\_{\overline{\pi}}} \exp \left[ -i \left( n\_{\mathcal{X}\overline{\mathcal{S}}} \overline{x} + n\_{yh} \overline{y} \right) \right] \\ &- \sum\_{g^h} -i n\_{\mathcal{X}\overline{\mathcal{S}}} e\_{z, \mathcal{S}h} \left( \overline{z}\_{\ell} \right) \exp \left[ -i \left( n\_{\mathcal{X}\overline{\mathcal{S}}} \overline{x} + n\_{yh} \overline{y} \right) \right] \\ &= -i \sqrt{Z\_0} H\_{\overline{\mathcal{Y}}} \\ &= -i \sum\_{g^h} h\_{\mathcal{Y}, \mathcal{S}h} \left( \overline{z}\_{\ell} \right) \exp \left[ -i \left( n\_{\mathcal{X}\overline{\mathcal{S}}} \overline{x} + n\_{yh} \overline{y} \right) \right] \end{split} \tag{59}$$

and is simplified to be:

$$\frac{\partial e\_{\rm x,gl}(\Xi\_{\ell})}{\partial \Xi} = -i\mathfrak{n}\_{\rm x\emptyset}e\_{\rm z,gl}(\Xi\_{\ell}) - i\mathfrak{h}\_{\rm y,gl}(\Xi\_{\ell}) \tag{60}$$

For the component *ez*,*gh*, the equation can be derived as:

*∂x <sup>Z</sup>*0*Hy* <sup>−</sup> *<sup>∂</sup><sup>y</sup> <sup>Z</sup>*0*Hx* <sup>=</sup> ∑ *gh* <sup>−</sup>*inxghy*,*gh* (*<sup>z</sup>*) exp −*i nxgx* + *nyhy* − ∑ *gh* <sup>−</sup>*inyhhx*,*gh* (*<sup>z</sup>*) exp −*i nxgx* + *nyhy* = *i <sup>Y</sup>*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**] *z* = *i* ∑ *ghuv <sup>ε</sup>zx*,*uvex*,*gh* (*<sup>z</sup>*) exp −*i nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y* +*i* ∑ *ghuv <sup>ε</sup>zy*,*uvey*,*gh* (*<sup>z</sup>*) exp −*i nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y* +*i* ∑ *ghuv <sup>ε</sup>zz*,*uvez*,*gh* (*<sup>z</sup>*) exp −*i nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y* (61)

and is simplified to be:

$$\begin{split} n\_{yh} h\_{\mathbf{x}, \operatorname{gl}} \left( \Xi\_{\ell} \right) - n\_{\mathbf{x}\mathcal{J}} h\_{\mathbf{y}, \operatorname{gl}} \left( \Xi\_{\ell} \right) &= \sum\_{\begin{subarray}{c} \mathsf{u}' \mathsf{v}' \\ \mathsf{u}' \mathsf{v}' \end{subarray}} \overline{\mathsf{c}}\_{\mathsf{zx}, (\mathsf{g} - \mathsf{u}') (\mathsf{h} - \mathsf{v}')} \mathsf{e}\_{\mathsf{x}, \mathsf{u}' \mathsf{v}'} \left( \Xi\_{\ell} \right) \\ &+ \sum\_{\begin{subarray}{c} \mathsf{u}' \mathsf{v}' \\ \mathsf{u}' \mathsf{v}' \end{subarray}} \overline{\mathsf{c}}\_{\mathsf{z}, (\mathsf{g} - \mathsf{u}') (\mathsf{h} - \mathsf{v}')} \mathsf{e}\_{\mathsf{y}, \mathsf{u}' \mathsf{v}'} \left( \Xi\_{\ell} \right) \\ &+ \sum\_{\begin{subarray}{c} \mathsf{u}' \mathsf{v}' \\ \mathsf{u}' \mathsf{v}' \end{subarray}} \overline{\mathsf{c}}\_{\mathsf{z}, \mathsf{u}' \mathsf{v}'} \left( \Xi\_{\ell} \right) \end{split} \tag{62}$$

For the component *<sup>∂</sup>hy*,*gh*(*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

$$\begin{split} \partial\_{\overline{y}}\sqrt{\varpi\_{0}}H\_{\overline{z}}-\partial\_{\overline{z}}\sqrt{\varpi\_{0}}H\_{y} &= \sum\_{g\mathcal{H}}-in\_{y\mathcal{H}}h\_{z,g\mathcal{H}}\left(\overline{\varpi}\_{\ell}\right)\exp\left[-i\left(n\_{xg}\overline{\varpi}+n\_{y\mathcal{H}}\overline{y}\right)\right] \\ &- \sum\_{m}\frac{\partial h\_{y,g\mathcal{H}}\left(\overline{\varpi}\_{\ell}\right)}{\partial\_{\overline{z}}}\exp\left[-i\left(n\_{xg}\overline{\varpi}+n\_{y\mathcal{H}}\overline{y}\right)\right] \\ &= i\sqrt{Y\_{0}}\left[\overline{\varpi}\left(\overline{r}\right)\to\right]\_{\mathtt{x}} \\ &= i\sum\_{g\mathcal{H}uv}\overline{\varepsilon}\_{\overline{x}x,uw}e\_{x,gh}\left(\overline{\varpi}\_{\ell}\right)\exp\left[-i\left(n\_{x\left(g+u\right)}\overline{\varpi}+n\_{y\left(h+v\right)}\overline{y}\right)\right] \\ &+ i\sum\_{g\mathcal{H}uv}\overline{\varepsilon}\_{x\mathcal{H},uw}e\_{z,gh}\left(\overline{\varpi}\_{\ell}\right)\exp\left[-i\left(n\_{x\left(g+u\right)}\overline{\varpi}+n\_{y\left(h+v\right)}\overline{y}\right)\right] \\ &+ i\sum\_{g\mathcal{H}uv}\overline{\varepsilon}\_{xx,uw}e\_{z,gh}\left(\overline{\varpi}\_{\ell}\right)\exp\left[-i\left(n\_{x\left(g+u\right)}\overline{\varpi}+n\_{y\left(h+v\right)}\overline{y}\right)\right] \end{split} \tag{63}$$

and is simplified to be:

12 Will-be-set-by-IN-TECH

*<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

*∂ex*,*gh* (*z*) *∂z*

*<sup>Z</sup>*0*Hy*

−*inxghy*,*gh* (*z*) exp

−*inyhhx*,*gh* (*z*) exp

*εzx*,*uvex*,*gh* (*z*) exp

*εzy*,*uvey*,*gh* (*z*) exp

*εzz*,*uvez*,*gh* (*z*) exp

*u*�*v*�

+ ∑ *u*�*v*�

+ ∑ *u*�*v*�

*<sup>ε</sup>zx*,(*g*−*u*�

*<sup>ε</sup>zy*,(*g*−*u*�

*<sup>ε</sup>zz*,(*g*−*u*�

exp −*i* 

−*inxgez*,*gh* (*z*) exp

*hy*,*gh* (*z*) exp

*m*

= −*i*

= −*i*∑ *gh*

− ∑ *gh*

= −*inyhez*,*gh* (*z*) + *ihx*,*gh* (*z*) (58)

*nxgx* + *nyhy*

*nxgx* + *nyhy*

= −*inxgez*,*gh* (*z*) − *ihy*,*gh* (*z*) (60)

*nxgx* + *nyhy*

*nxgx* + *nyhy*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*nx*(*g*+*u*)*x* + *ny*(*h*+*v*)*y*

*<sup>v</sup>*� (*z*)

*<sup>v</sup>*� (*z*)

 −*i* 

 −*i* 

 −*i* 

> −*i*

 −*i* 

> −*i*

 −*i* 

)(*h*−*v*�

)(*h*−*v*�

)(*h*−*v*�

)*ex*,*u*�

)*ey*,*u*�

)*ez*,*u*�

*nxgx* + *nyhy*

*<sup>v</sup>*� (*z*) (62)

(61)

(59)

*∂ey*,*gh* (*z*) *∂z*

*<sup>Y</sup>*0*Ez* <sup>=</sup> ∑

*∂ex*,*gh* (*z*) *∂z*

For the component *ez*,*gh*, the equation can be derived as:

*<sup>Z</sup>*0*Hx* <sup>=</sup> ∑

*gh*

= *i*

= *i* ∑ *ghuv*

> +*i* ∑ *ghuv*

> +*i* ∑ *ghuv*

*nyhhx*,*gh* (*<sup>z</sup>*) − *nxghy*,*gh* (*<sup>z</sup>*) = ∑

− ∑ *gh*

*<sup>Y</sup>*<sup>0</sup> [*ε*(*r*) **<sup>E</sup>**]*<sup>z</sup>*

and is simplified to be:

For the component *<sup>∂</sup>ex*,*gh* (*<sup>z</sup>* )

*∂z*

and is simplified to be:

*∂x*

and is simplified to be:

*<sup>Z</sup>*0*Hy* <sup>−</sup> *<sup>∂</sup><sup>y</sup>*

*<sup>Y</sup>*0*Ex* <sup>−</sup> *<sup>∂</sup><sup>x</sup>*

$$\begin{split} \frac{\partial h\_{y,\text{gf}}(\overline{\boldsymbol{\tau}}\_{\ell})}{\partial \overline{\boldsymbol{\tau}}} &= -i\mathfrak{n}\_{y\boldsymbol{h}}h\_{\overline{\boldsymbol{z}},\text{gf}}\left(\overline{\boldsymbol{\tau}}\_{\ell}\right) - i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\varepsilon}}\_{\text{xx},(\underline{\boldsymbol{g}}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\mathfrak{e}\_{\boldsymbol{x},\boldsymbol{u}'\boldsymbol{v}'}\left(\overline{\boldsymbol{\varpi}}\_{\ell}\right) \\ &-i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\varepsilon}}\_{\text{xy},(\underline{\boldsymbol{g}}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\mathfrak{e}\_{\boldsymbol{y},\boldsymbol{u}'\boldsymbol{v}'}\left(\overline{\boldsymbol{\varpi}}\_{\ell}\right) - i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\varepsilon}}\_{\text{xz},(\underline{\boldsymbol{g}}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\mathfrak{e}\_{\boldsymbol{z},\boldsymbol{u}'\boldsymbol{v}'}\left(\overline{\boldsymbol{\varpi}}\_{\ell}\right) \end{split} \tag{64}$$

For the component *<sup>∂</sup>hx*,*gh* (*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , the equation can be derived as:

$$\begin{split} \partial\_{\overline{\Xi}} \sqrt{Z\_{0}} H\_{\text{X}} - \partial\_{\overline{\Xi}} \sqrt{Z\_{0}} H\_{\text{Z}} &= \sum\_{gh} \frac{\partial t\_{x,gh} \left( \overline{z\_{\ell}} \right)}{\partial\_{\overline{\Xi}}} \exp \left[ -i \left( n\_{x\overline{\xi}} \overline{x} + n\_{yh} \overline{y} \right) \right] \\ &- \sum\_{gh} -i n\_{xgh} n\_{z,gh} \left( \overline{z\_{\ell}} \right) \exp \left[ -i \left( n\_{x\overline{\xi}} \overline{x} + n\_{yh} \overline{y} \right) \right] \\ &= i \sqrt{Y\_{0}} \left[ \overline{e} \left( \overline{\tau} \right) \mathbf{E} \right]\_{g} \\ &= i \sum\_{g \text{huv}} \overline{e}\_{yx, \mu v} e\_{x, gh} \left( \overline{\tau\_{\ell}} \right) \exp \left[ -i \left( n\_{x(\underline{\xi} + \underline{u})} \overline{\mathbf{x}} + n\_{y(h + v)} \overline{y} \right) \right] \\ &+ i \sum\_{g \text{huv}} \overline{e}\_{yy, \mu v} e\_{z, gh} \left( \overline{\tau\_{\ell}} \right) \exp \left[ -i \left( n\_{x(\underline{\xi} + \underline{u})} \overline{\mathbf{x}} + n\_{y(h + v)} \overline{y} \right) \right] \\ &+ i \sum\_{g \text{huv}} \overline{e}\_{yz, \mu v} e\_{z, gh} \left( \overline{\tau\_{\ell}} \right) \exp \left[ -i \left( n\_{x(\underline{\xi} + \underline{u})} \overline{\mathbf{x}} + n\_{y(h + v)} \overline{y} \right) \right] \tag{65} \end{split} \tag{66}$$

and is simplified to be:

$$\begin{split} \frac{\partial h\_{\mathbf{x},\mathbf{g},\mathbf{h}}(\Xi\_{\ell})}{\partial \Xi} &= -i\mathbf{n}\_{\mathbf{x}\mathbf{h}}h\_{\mathbf{z},\mathbf{g},\mathbf{h}}(\Xi\_{\ell}) + i\sum\_{\mathbf{u}'\mathbf{v}'}\mathsf{\mathbf{\tilde{z}}}\_{\mathbf{y}\mathbf{x},(\mathbf{g}-\mathbf{u}')(\mathbf{h}-\mathbf{v}')}\mathfrak{e}\_{\mathbf{x},\mathbf{u}'\mathbf{v}'}(\Xi\_{\ell}) \\ &+ i\sum\_{\mathbf{u}'\mathbf{v}'}\mathsf{\tilde{z}}\_{\mathbf{y}\mathbf{y},(\mathbf{g}-\mathbf{u}')(\mathbf{h}-\mathbf{v}')}\mathfrak{e}\_{\mathbf{y},\mathbf{u}'\mathbf{v}'}(\Xi\_{\ell}) + i\sum\_{\mathbf{u}'\mathbf{v}'}\mathsf{\tilde{z}}\_{\mathbf{y}\mathbf{z},(\mathbf{g}-\mathbf{u}')(\mathbf{h}-\mathbf{v}')}\mathfrak{e}\_{\mathbf{z},\mathbf{u}'\mathbf{v}'}(\Xi\_{\ell}) \end{split} \tag{66}$$

Here, the symbol (.) represents a *NgNh* <sup>×</sup> 1 vector, and the symbol ˜(.) represents a *NgNh* <sup>×</sup> *NgNh* matrix, indicating the considered *g*(*h*) ranged from *gmin*(*hmin*) to *gmax*(*hmax*) with *Ng* = *gmin* + *gmax* + 1(*Nh* = *hmin* + *hmax* + 1). As indicated in Equation (69), the normal field **f***n*<sup>ˆ</sup> can

<sup>77</sup> Electromagnetic Formalisms for Optical Propagation

Further, we derive the coupled-wave equation of the tangential field **f**ˆ*t*, in which the component fields of **f***n*<sup>ˆ</sup> are replaced by those of **f**ˆ*<sup>t</sup>* via equation (69). Similarly, we consider

= −*inxgez*,*gh* (*z*) − *ihy*,*gh* (*z*) (60)

= −*inyhez*,*gh* (*z*) + *ihx*,*gh* (*z*) (58)

)(*h*−*v*�

)(*h*−*v*�

*u*� *v*�

*u*�*v*�

*<sup>v</sup>*� (*<sup>z</sup>*) − *<sup>i</sup>* ∑

*<sup>v</sup>*� (*z*) + *i* ∑

−*n*˜ *<sup>x</sup>ez*

*hz* <sup>−</sup> *<sup>ε</sup>*˜*xzez* −*n*˜ *<sup>y</sup>ez*

*hz* + *ε*˜*yzez*

−1

−1

*zz n*˜ *<sup>x</sup>* − 1 *n*˜ *<sup>x</sup>ε*˜

*zz n*˜ *<sup>x</sup> n*˜ *<sup>y</sup>ε*˜

·*i***f**ˆ*<sup>t</sup>* ≡ *i***G** · **f**ˆ*<sup>t</sup>* (70)

*zz n*˜ *<sup>x</sup> ε*˜*xzε*˜

*zz n*˜ *<sup>x</sup>* −*ε*˜*yzε*˜

Definitely, the equation (70) turns the Maxwell's curl equations into a eigen-system problems. Up to now, with the known structured layers for equations (44)-(45) and the known incidence related to equations (46)-(47), the transition behaviors of the tangential field **f**ˆ*<sup>t</sup>* can be formulated layer by layer via equation (70), and the corresponding normal field **f***n*<sup>ˆ</sup> can be

*<sup>∂</sup><sup>z</sup>* , and *<sup>∂</sup>hx*,*gh* (*<sup>z</sup>* )

)*ex*,*u*�

)*ex*,*u*�

*<sup>v</sup>*� (*z*)

*<sup>v</sup>*� (*z*)

)(*h*−*v*�

)(*h*−*v*�

)*ez*,*u*�

)*ez*,*u*�

*zz ε*˜*zy* −*n*˜ *<sup>x</sup>ε*˜

*zz ε*˜*zy* −*n*˜ *<sup>y</sup>ε*˜

*zz <sup>ε</sup>*˜*zy* <sup>+</sup> *<sup>ε</sup>*˜*yy* <sup>−</sup> *<sup>n</sup>*˜ *xn*˜ *<sup>x</sup> <sup>ε</sup>*˜*yzε*˜−<sup>1</sup>

*zz ε*˜*zy* − *ε*˜*xy* − *n*˜ *yn*˜ *<sup>x</sup>* −*ε*˜*xzε*˜

*<sup>v</sup>*� (*z*) (64)

*<sup>v</sup>*� (*z*) (66)

−1 *zz n*˜ *<sup>y</sup>* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

−1 *zz n*˜ *<sup>y</sup>*

*zz n*˜ *<sup>y</sup>*

−1 *zz n*˜ *<sup>y</sup>* + 1

*<sup>ε</sup>xz*,(*g*−*u*�

*<sup>ε</sup>yz*,(*g*−*u*�

⎤ ⎥ ⎥ ⎥ ⎦

−1

−1

*<sup>∂</sup><sup>z</sup>* in Equations (58), (60), (64),

*<sup>∂</sup><sup>z</sup>* , *<sup>∂</sup>hy*,*gh* (*<sup>z</sup>* )

*<sup>ε</sup>xx*,(*g*−*u*�

*<sup>ε</sup>yx*,(*g*−*u*�

be obtained straightforwardly if the tangential field **f**ˆ*<sup>t</sup>* is given.

in Three-Dimensional Periodic Liquid-Crystal Microstructures

**A.4 Derive the coupled-wave equation of the tangential field f**ˆ*<sup>t</sup>*

= −*inyhhz*,*gh* (*<sup>z</sup>*) − *<sup>i</sup>* ∑

*<sup>ε</sup>xy*,(*g*−*u*�

= −*inxhhz*,*gh* (*<sup>z</sup>*) + *<sup>i</sup>* ∑

*<sup>ε</sup>yy*,(*g*−*u*�

⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*zz ε*˜*zx n*˜ *<sup>x</sup>ε*˜

*zz ε*˜*zx n*˜ *<sup>y</sup>ε*˜

<sup>−</sup><sup>1</sup> *zz <sup>ε</sup>*˜*zx* <sup>−</sup> *<sup>ε</sup>*˜*xx* <sup>+</sup> *<sup>n</sup>*˜ *yn*˜ *<sup>y</sup> <sup>ε</sup>*˜*xzε*˜

*zz ε*˜*zx* + *ε*˜*yx* + *n*˜ *xn*˜ *<sup>y</sup>* −*ε*˜*yzε*˜

*ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + *i*

−*<sup>i</sup>* ∑ *u*�*v*�

+*i* ∑ *u*�*v*�

0 −100 −*ε*˜*xx* 0 −*ε*˜*xy* 0 0 0 01 *ε*˜*yx* 0 *ε*˜*yy* 0

> *n*˜ *<sup>x</sup>ε*˜ −1

> *n*˜ *yε*˜ −1

evaluated sequentially via equation (69).

*<sup>∂</sup><sup>z</sup>* , *<sup>∂</sup>ex*,*gh* (*<sup>z</sup>* )

*u*�*v*�

*u*� *v*�

)*ey*,*u*�

)*ey*,*u*�

With equation (69), these equations can matrixize the coupled-wave equation of **f**ˆ*<sup>t</sup>* as:

⎡ ⎢ ⎢ ⎢ ⎣

−1

−1

−1

−1

−*n*˜ *<sup>y</sup>*

−*n*˜ *<sup>x</sup>*

)(*h*−*v*�

)(*h*−*v*�

the associated formulas of *<sup>∂</sup>ey*,*gh*(*<sup>z</sup>* )

and (66), respectively, i.e.:

*∂ex*,*gh* (*z*) *∂z*

*∂ey*,*gh* (*z*) *∂z*

*∂hy*,*gh* (*z*) *∂z*

*∂hx*,*gh* (*z*) *∂z*

> ⎡ ⎢ ⎢ ⎣

*∂***f**ˆ*t ∂z* = *i*

=

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*ε*˜*xzε*˜

<sup>−</sup>*ε*˜*yzε*˜−<sup>1</sup>

To solve the fields systematically, these equations are reformulated in terms of the full fields **f**ˆ*<sup>t</sup>* and **f***n*<sup>ˆ</sup> in the following context, and show an eigen-system problem for studied periodic structures.

#### **A.3 Derive the coupled-wave equation of the normal field f***n*<sup>ˆ</sup>

To obtain the coupled-wave equations of the normal field **f***n*ˆ, we consider the above-mentioned formulas for its components *hz*,*gh* (*z*) and *ez*,*gh* (*z*) in Equations (56) and (62), respectively, i.e.:

$$h\_{z,gh}(\overline{\varpi}\_{\ell}) = n\_{\overline{x}\overline{g}}e\_{y,gh}(\overline{\varpi}\_{\ell}) - n\_{yh}e\_{x,gh}(\overline{\varpi}\_{\ell}) \tag{56}$$

$$\begin{split} n\_{yh}h\_{x,gh}(\overline{\varpi}\_{\ell}) - n\_{\overline{x}\overline{g}}h\_{y,gh}(\overline{\varpi}\_{\ell}) &= \sum\_{\mathsf{u}'\mathsf{v}'} \overline{\mathsf{E}}\_{\mathsf{zx},(\underline{g}-\mathsf{u}')(h-\mathsf{v}')} e\_{\mathsf{x},\mathsf{u}'\mathsf{v}'}(\overline{\varpi}\_{\ell}) \\ &+ \sum\_{\mathsf{u}'\mathsf{v}'} \overline{\mathsf{E}}\_{\overline{x}\underline{y},(\underline{g}-\mathsf{u}')(h-\mathsf{v}')} e\_{\mathsf{y},\mathsf{u}'\mathsf{v}'}(\overline{\varpi}\_{\ell}) \\ &+ \sum\_{\mathsf{u}'\mathsf{v}'} \overline{\mathsf{E}}\_{\overline{x}\underline{z},(\underline{g}-\mathsf{u}')(h-\mathsf{v}')} e\_{\mathsf{z},\mathsf{u}'\mathsf{v}'}(\overline{\varpi}\_{\ell}) \end{split} \tag{62}$$

Up to the Fourier order *g*, *h* ∈ {0, 1}, an example corresponding to Equations (56) and (62) can be matrixized:

⎡ ⎢ ⎢ ⎢ ⎣ *hz*,00 (*z*) *hz*,01 (*z*) *hz*,10 (*z*) *hz*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ *nx*<sup>0</sup> 000 0 *nx*<sup>0</sup> 0 0 0 0 *nx*<sup>1</sup> 0 000 *nx*<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ey*,00 (*z*) *ey*,01 (*z*) *ey*,10 (*z*) *ey*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ − ⎡ ⎢ ⎢ ⎢ ⎣ *ny*<sup>0</sup> 000 0 *ny*<sup>1</sup> 0 0 0 0 *ny*<sup>0</sup> 0 000 *ny*<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ex*,00 (*z*) *ex*,01 (*z*) *ex*,10 (*z*) *ex*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ (67) ⎡ ⎢ ⎢ ⎢ ⎣ *εzz*,00 *εzz*,0−<sup>1</sup> *εzz*,−<sup>10</sup> *εzz*,−1−<sup>1</sup> *εzz*,01 *εzz*,00 *εzz*,−<sup>11</sup> *εzz*,−<sup>10</sup> *εzz*,10 *εzz*,1−<sup>1</sup> *εzz*,00 *εzz*,0−<sup>1</sup> *εzz*,11 *εzz*,10 *εzz*,01 *εzz*,00 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ez*,00 (*z*) *ez*,01 (*z*) *ez*,10 (*z*) *ez*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ *ny*<sup>0</sup> 000 0 *ny*<sup>1</sup> 0 0 0 0 *ny*<sup>0</sup> 0 000 *ny*<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *hx*,00 (*z*) *hx*,01 (*z*) *hx*,10 (*z*) *hx*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ − ⎡ ⎢ ⎢ ⎢ ⎣ *nx*<sup>0</sup> 000 0 *nx*<sup>0</sup> 0 0 0 0 *nx*<sup>1</sup> 0 000 *nx*<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *hy*,00 (*z*) *hy*,01 (*z*) *hy*,10 (*z*) *hy*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ − ⎡ ⎢ ⎢ ⎢ ⎣ *εzx*,00 *εzx*,0−<sup>1</sup> *εzx*,−<sup>10</sup> *εzx*,−1−<sup>1</sup> *εzx*,01 *εzx*,00 *εzx*,−<sup>11</sup> *εzx*,−<sup>10</sup> *εzx*,10 *εzx*,1−<sup>1</sup> *εzx*,00 *εzx*,0−<sup>1</sup> *εzx*,11 *εzx*,10 *εzx*,01 *εzx*,00 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ex*,00 (*z*) *ex*,01 (*z*) *ex*,10 (*z*) *ex*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ − ⎡ ⎢ ⎢ ⎢ ⎣ *εzy*,00 *εzy*,0−<sup>1</sup> *εzy*,−<sup>10</sup> *εzy*,−1−<sup>1</sup> *εzy*,01 *εzy*,00 *εzy*,−<sup>11</sup> *εzy*,−<sup>10</sup> *εzy*,10 *εzy*,1−<sup>1</sup> *εzy*,00 *εzy*,0−<sup>1</sup> *εzy*,11 *εzy*,10 *εzy*,01 *εzy*,00 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ey*,00 (*z*) *ey*,01 (*z*) *ey*,10 (*z*) *ey*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ (68)

The full-component coupled-wave equation for the normal field **f***n*<sup>ˆ</sup> then can be extended as:

$$\mathbf{f}\_{\boldsymbol{\theta}} = \begin{bmatrix} \vec{e}\_{z} \\ \vec{h}\_{z} \end{bmatrix} = \begin{bmatrix} -\vec{\varepsilon}\_{zz}^{-1}\vec{\varepsilon}\_{xx} - \vec{\varepsilon}\_{zz}^{-1}\vec{n}\_{x} - \vec{\varepsilon}\_{zz}^{-1}\vec{\varepsilon}\_{zy}\,\,\vec{\varepsilon}\_{zz}^{-1}\vec{n}\_{y} \\\ -\vec{n}\_{y} & 0 \qquad \vec{n}\_{x} & 0 \end{bmatrix} \cdot \begin{bmatrix} \vec{e}\_{x} \\ \vec{h}\_{y} \\ \vec{e}\_{y} \\ \vec{h}\_{x} \end{bmatrix} \equiv \mathbf{D} \cdot \mathbf{f}\_{\hat{l}} \tag{69}$$

14 Will-be-set-by-IN-TECH

To solve the fields systematically, these equations are reformulated in terms of the full fields **f**ˆ*<sup>t</sup>* and **f***n*<sup>ˆ</sup> in the following context, and show an eigen-system problem for studied periodic

To obtain the coupled-wave equations of the normal field **f***n*ˆ, we consider the above-mentioned formulas for its components *hz*,*gh* (*z*) and *ez*,*gh* (*z*) in Equations (56) and

*<sup>ε</sup>zx*,(*g*−*u*�

*<sup>ε</sup>zy*,(*g*−*u*�

*<sup>ε</sup>zz*,(*g*−*u*�

Up to the Fourier order *g*, *h* ∈ {0, 1}, an example corresponding to Equations (56) and (62) can

*ey*,00 (*z*) *ey*,01 (*z*) *ey*,10 (*z*) *ey*,11 (*z*)

*ez*,00 (*z*) *ez*,01 (*z*) *ez*,10 (*z*) *ez*,11 (*z*)

> ⎡ ⎢ ⎢ ⎢ ⎣

*hz*,*gh* (*z*) = *nxgey*,*gh* (*z*) − *nyhex*,*gh* (*z*) (56)

)*ex*,*u*�

)*ey*,*u*�

)*ez*,*u*�

*<sup>v</sup>*� (*z*)

*<sup>v</sup>*� (*z*)

*ny*<sup>0</sup> 000 0 *ny*<sup>1</sup> 0 0 0 0 *ny*<sup>0</sup> 0 000 *ny*<sup>1</sup>

*ny*<sup>0</sup> 000 0 *ny*<sup>1</sup> 0 0 0 0 *ny*<sup>0</sup> 0 000 *ny*<sup>1</sup>

*εzx*,00 *εzx*,0−<sup>1</sup> *εzx*,−<sup>10</sup> *εzx*,−1−<sup>1</sup> *εzx*,01 *εzx*,00 *εzx*,−<sup>11</sup> *εzx*,−<sup>10</sup> *εzx*,10 *εzx*,1−<sup>1</sup> *εzx*,00 *εzx*,0−<sup>1</sup> *εzx*,11 *εzx*,10 *εzx*,01 *εzx*,00

> � ·

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

*<sup>v</sup>*� (*z*) (62)

⎤ ⎥ ⎥ ⎥ ⎦

⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣

> ⎤ ⎥ ⎥ ⎥ ⎦

⎡ ⎢ ⎢ ⎢ ⎣ *ex*,00 (*z*) *ex*,01 (*z*) *ex*,10 (*z*) *ex*,11 (*z*)

≡ **D** · **f**ˆ*<sup>t</sup>* (69)

⎤ ⎥ ⎥ ⎥ ⎦

(68)

⎡ ⎢ ⎢ ⎢ ⎣ *ex*,00 (*z*) *ex*,01 (*z*) *ex*,10 (*z*) *ex*,11 (*z*)

*hx*,00 (*z*) *hx*,01 (*z*) *hx*,10 (*z*) *hx*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦

⎤ ⎥ ⎥ ⎥ ⎦ (67)

)(*h*−*v*�

)(*h*−*v*�

)(*h*−*v*�

⎤ ⎥ ⎥ ⎥ ⎦ −

⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎣

**A.3 Derive the coupled-wave equation of the normal field f***n*<sup>ˆ</sup>

*u*� *v*�

+ ∑ *u*� *v*�

+ ∑ *u*�*v*�

⎤ ⎥ ⎥ ⎥ ⎦

⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣

> ⎤ ⎥ ⎥ ⎥ ⎦ −

⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ *ey*,00 (*z*) *ey*,01 (*z*) *ey*,10 (*z*) *ey*,11 (*z*)

The full-component coupled-wave equation for the normal field **f***n*<sup>ˆ</sup> then can be extended as:

−*n*˜ *<sup>y</sup>* 0 *n*˜ *<sup>x</sup>* 0

−1 *zz n*˜ *<sup>x</sup>* −*ε*˜ ⎤ ⎥ ⎥ ⎥ ⎦

−1 *zz ε*˜*zy ε*˜ −1 *zz n*˜ *<sup>y</sup>*

⎡ ⎢ ⎢ ⎢ ⎣

*nyhhx*,*gh* (*<sup>z</sup>*) − *nxghy*,*gh* (*<sup>z</sup>*) = ∑

structures.

(62), respectively, i.e.:

be matrixized: ⎡ ⎢ ⎢ ⎢ ⎣

> ⎡ ⎢ ⎢ ⎢ ⎣

−

−

⎡ ⎢ ⎢ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎣

*hz*,00 (*z*) *hz*,01 (*z*) *hz*,10 (*z*) *hz*,11 (*z*) ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣

*εzz*,00 *εzz*,0−<sup>1</sup> *εzz*,−<sup>10</sup> *εzz*,−1−<sup>1</sup> *εzz*,01 *εzz*,00 *εzz*,−<sup>11</sup> *εzz*,−<sup>10</sup> *εzz*,10 *εzz*,1−<sup>1</sup> *εzz*,00 *εzz*,0−<sup>1</sup> *εzz*,11 *εzz*,10 *εzz*,01 *εzz*,00

> ⎤ ⎥ ⎥ ⎥ ⎦

*εzy*,00 *εzy*,0−<sup>1</sup> *εzy*,−<sup>10</sup> *εzy*,−1−<sup>1</sup> *εzy*,01 *εzy*,00 *εzy*,−<sup>11</sup> *εzy*,−<sup>10</sup> *εzy*,10 *εzy*,1−<sup>1</sup> *εzy*,00 *εzy*,0−<sup>1</sup> *εzy*,11 *εzy*,10 *εzy*,01 *εzy*,00

⎡ ⎢ ⎢ ⎢ ⎣ *hy*,00 (*z*) *hy*,01 (*z*) *hy*,10 (*z*) *hy*,11 (*z*)

*nx*<sup>0</sup> 000 0 *nx*<sup>0</sup> 0 0 0 0 *nx*<sup>1</sup> 0 000 *nx*<sup>1</sup>

> **f***n*<sup>ˆ</sup> = � *ez hz* � = � −*ε*˜ −1 *zz ε*˜*zx* −*ε*˜

*nx*<sup>0</sup> 000 0 *nx*<sup>0</sup> 0 0 0 0 *nx*<sup>1</sup> 0 000 *nx*<sup>1</sup> Here, the symbol (.) represents a *NgNh* <sup>×</sup> 1 vector, and the symbol ˜(.) represents a *NgNh* <sup>×</sup> *NgNh* matrix, indicating the considered *g*(*h*) ranged from *gmin*(*hmin*) to *gmax*(*hmax*) with *Ng* = *gmin* + *gmax* + 1(*Nh* = *hmin* + *hmax* + 1). As indicated in Equation (69), the normal field **f***n*<sup>ˆ</sup> can be obtained straightforwardly if the tangential field **f**ˆ*<sup>t</sup>* is given.

## **A.4 Derive the coupled-wave equation of the tangential field f**ˆ*<sup>t</sup>*

Further, we derive the coupled-wave equation of the tangential field **f**ˆ*t*, in which the component fields of **f***n*<sup>ˆ</sup> are replaced by those of **f**ˆ*<sup>t</sup>* via equation (69). Similarly, we consider the associated formulas of *<sup>∂</sup>ey*,*gh*(*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , *<sup>∂</sup>ex*,*gh* (*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , *<sup>∂</sup>hy*,*gh* (*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* , and *<sup>∂</sup>hx*,*gh* (*<sup>z</sup>* ) *<sup>∂</sup><sup>z</sup>* in Equations (58), (60), (64), and (66), respectively, i.e.:

$$\frac{\partial e\_{x,glt}(\overline{z}\_{\ell})}{\partial \overline{z}\_{\overline{z}}} = -i\imath\_{X\overline{g}} e\_{z,glt}(\overline{z}\_{\ell}) - i\hbar\_{\overline{g},glt}(\overline{z}\_{\ell}) \tag{60}$$

$$\frac{\partial \varepsilon\_{y,gl}(\Xi\_{\ell})}{\partial \Xi} = -i n\_{yh} \varepsilon\_{z,gl}(\Xi\_{\ell}) + i h\_{x,gl}(\Xi\_{\ell}) \tag{58}$$

$$\begin{split} \frac{\partial h\_{y,\text{ghl}}(\overline{\boldsymbol{\tau}}\_{\ell})}{\partial \overline{\boldsymbol{\tau}}\_{\ell}} &= -i\boldsymbol{n}\_{y\boldsymbol{h}}h\_{z,\text{gph}}(\overline{\boldsymbol{\tau}}\_{\ell}) - i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\tau}}\_{\text{xx},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{e}\_{\text{x},\text{h}'\text{v}'}\left(\overline{\boldsymbol{\tau}}\_{\ell}\right) \\ &-i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\tau}}\_{\text{xy},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{e}\_{\text{y},\text{h}'\text{v}'}\left(\overline{\boldsymbol{\tau}}\_{\ell}\right) - i\sum\_{\boldsymbol{u}'\boldsymbol{v}'}\overline{\boldsymbol{\tau}}\_{\text{xz},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{e}\_{\text{z},\text{h}'\text{v}'}\left(\overline{\boldsymbol{\tau}}\_{\ell}\right) \\ \partial\boldsymbol{h}\_{\text{x},\text{gph}}\left(\overline{\boldsymbol{\tau}}\_{\ell}\right) &=&\operatorname{i}\boldsymbol{u}\cdot\boldsymbol{k} \qquad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \quad \boldsymbol{\tau}\_{\ell} \end{split} \tag{64}$$

$$\begin{split} \frac{\partial \mathbf{u}\_{\mathbf{x},\mathbf{g}h}(\boldsymbol{\omega}\_{\ell})}{\partial \mathbf{z}} &= -i\mathbf{n}\_{\mathbf{x}}h\_{\mathbf{z},\mathbf{g}h}(\boldsymbol{\Xi}\_{\ell}) + i\sum\_{\mathbf{u}'\mathbf{v}'}\boldsymbol{\overline{\varepsilon}}\_{\mathbf{y}\mathbf{x},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{\varepsilon}\_{\mathbf{x},\mathbf{u}'\mathbf{v}'}(\boldsymbol{\Xi}\_{\ell}) \\ &+ i\sum\_{\mathbf{u}'\mathbf{v}'}\boldsymbol{\overline{\varepsilon}}\_{\mathbf{y}\mathbf{y},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{\varepsilon}\_{\mathbf{y},\mathbf{u}'\mathbf{v}'}(\boldsymbol{\Xi}\_{\ell}) + i\sum\_{\mathbf{u}'\mathbf{v}'}\boldsymbol{\overline{\varepsilon}}\_{\mathbf{y}\mathbf{z},(\boldsymbol{g}-\boldsymbol{u}')(\boldsymbol{h}-\boldsymbol{v}')}\boldsymbol{\varepsilon}\_{\mathbf{z},\mathbf{u}'\mathbf{v}'}(\boldsymbol{\Xi}\_{\ell}) \end{split} \tag{66}$$

With equation (69), these equations can matrixize the coupled-wave equation of **f**ˆ*<sup>t</sup>* as:

*∂***f**ˆ*t ∂z* = *i* ⎡ ⎢ ⎢ ⎣ 0 −100 −*ε*˜*xx* 0 −*ε*˜*xy* 0 0 0 01 *ε*˜*yx* 0 *ε*˜*yy* 0 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + *i* ⎡ ⎢ ⎢ ⎢ ⎣ −*n*˜ *<sup>x</sup>ez* −*n*˜ *<sup>y</sup> hz* <sup>−</sup> *<sup>ε</sup>*˜*xzez* −*n*˜ *<sup>y</sup>ez* −*n*˜ *<sup>x</sup> hz* + *ε*˜*yzez* ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *n*˜ *<sup>x</sup>ε*˜ −1 *zz ε*˜*zx n*˜ *<sup>x</sup>ε*˜ −1 *zz n*˜ *<sup>x</sup>* − 1 *n*˜ *<sup>x</sup>ε*˜ −1 *zz ε*˜*zy* −*n*˜ *<sup>x</sup>ε*˜ −1 *zz n*˜ *<sup>y</sup> ε*˜*xzε*˜ −1 *zz ε*˜*zx* − *ε*˜*xx* + *n*˜ *yn*˜ *<sup>y</sup> ε*˜*xzε*˜ −1 *zz n*˜ *<sup>x</sup> ε*˜*xzε*˜−<sup>1</sup> *zz <sup>ε</sup>*˜*zy* <sup>−</sup> *<sup>ε</sup>*˜*xy* <sup>−</sup> *<sup>n</sup>*˜ *yn*˜ *<sup>x</sup>* <sup>−</sup>*ε*˜*xzε*˜−<sup>1</sup> *zz n*˜ *<sup>y</sup> n*˜ *yε*˜ −1 *zz ε*˜*zx n*˜ *<sup>y</sup>ε*˜ −1 *zz n*˜ *<sup>x</sup> n*˜ *<sup>y</sup>ε*˜ −1 *zz ε*˜*zy* −*n*˜ *<sup>y</sup>ε*˜ −1 *zz n*˜ *<sup>y</sup>* + 1 <sup>−</sup>*ε*˜*yzε*˜−<sup>1</sup> *zz ε*˜*zx* + *ε*˜*yx* + *n*˜ *xn*˜ *<sup>y</sup>* −*ε*˜*yzε*˜ −1 *zz n*˜ *<sup>x</sup>* −*ε*˜*yzε*˜ −1 *zz ε*˜*zy* + *ε*˜*yy* − *n*˜ *xn*˜ *<sup>x</sup> ε*˜*yzε*˜ −1 *zz n*˜ *<sup>y</sup>* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ·*i***f**ˆ*<sup>t</sup>* ≡ *i***G** · **f**ˆ*<sup>t</sup>* (70)

Definitely, the equation (70) turns the Maxwell's curl equations into a eigen-system problems. Up to now, with the known structured layers for equations (44)-(45) and the known incidence related to equations (46)-(47), the transition behaviors of the tangential field **f**ˆ*<sup>t</sup>* can be formulated layer by layer via equation (70), and the corresponding normal field **f***n*<sup>ˆ</sup> can be evaluated sequentially via equation (69).

with *<sup>ξ</sup>gh* <sup>=</sup> �*<sup>ε</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg* while the corresponding eigen-vector matrix are:

*gh*<sup>3</sup> **v**� *gh*4 �

*nxgnxg*−*ε gh nxgnyh*

<sup>79</sup> Electromagnetic Formalisms for Optical Propagation

−*ε ghξ gh nxgnyh*

*ξ gh nxgnyh*

*nyhnyh*−*ε gh nxgnyh*

0101 1010

*gh*<sup>1</sup> − *nxg***v**�

*gh*<sup>1</sup> + *nyh***v**�

*gh*<sup>3</sup> − *nxg***v**�

*gh*<sup>3</sup> + *nyh***v**�

**v***gh*<sup>1</sup> **v***gh*<sup>2</sup> **v***gh*<sup>3</sup> **v***gh*<sup>4</sup>

**v**�

*gh*2 �

> *gh*2 �

*gh*4 �

*gh*4 �

*nyh mgh*

−*nyhξ gh mgh*

−*nxg mgh*

−*nxg ξ gh mgh*

�

*nxg mgh* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

ˆ*j* (80)

*k* (81)

<sup>−</sup>*<sup>ε</sup> ghnxg <sup>ξ</sup>*−<sup>1</sup> *gh mgh*

*nyh mgh*

*ε ghnyhξ*−<sup>1</sup> *gh mgh*

*nxgnxg*−*ε gh nxgnyh*

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

*gh*1, **v**�

/*mgh* (75)

/*mgh* (76)

/*mgh* (77)

/*mgh* (78)

(73)

(79)

*gh*2) as well as

*ε gh ξ gh nxgnyh*

*nyhnyh* + *nxgnxg* (74)

*gh*<sup>2</sup> **v**�

−*ξ gh nxgnyh*

*nyhnyh*−*ε gh nxgnyh*

Due to the degeneracy in (*κgh*,1, *κgh*,2) and (*κgh*,3, *κgh*,4), the eigenvector (**v**�

*gh*4) can be remixed by arbitrary linear combinations. Choosing

*mgh* <sup>=</sup> �

*nxgξgh***v**�

<sup>−</sup> *<sup>ε</sup>ghnyh ξgh*

−*nxgξgh***v**�

**v**�

*nxg mgh*

*ε ghnxg ξ*−<sup>1</sup> *gh mgh*

*nyh mgh*

<sup>−</sup>*<sup>ε</sup> ghnyhξ*−<sup>1</sup> *gh mgh*

Hence, **v***gh*<sup>1</sup> and **v***gh*<sup>2</sup> correspond to the (*g*, *h*)-order forward TE and TM (transverse electric and transverse magnetic) representations (with respect to the plane of the diffraction wave), respectively. **v***gh*<sup>3</sup> and **v***gh*<sup>4</sup> then correspond to the backward TE and TM ones. For example,

> *<sup>ı</sup>*ˆ<sup>−</sup> *nxg mgh*

> > *ı*ˆ+

*nyhξgh mgh*

<sup>ˆ</sup>*<sup>j</sup>* <sup>−</sup> *mgh* <sup>ˆ</sup>

� *<sup>ε</sup>ghnyh ξgh*

*eigvec* =

(**v**� *gh*3, **v**� � **v**� *gh*<sup>1</sup> **v**�

in Three-Dimensional Periodic Liquid-Crystal Microstructures

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

**v***gh*<sup>1</sup> =

**v***gh*<sup>2</sup> =

**v***gh*<sup>3</sup> =

**v***gh*<sup>4</sup> =

*eigvec* = **<sup>T</sup>***gh* = �

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

*nyh mgh*

*nyhξ gh mgh*

−*nxg mgh*

*nxg ξ gh mgh*

=

with equation (69) and (79), **v***gh*<sup>1</sup> denotes the component fields:

**<sup>e</sup>***gh* <sup>=</sup> *nyh mgh*

**<sup>h</sup>***gh* <sup>=</sup> *nxgξgh mgh*

the equation (73) is then shown as:

�

�

�

=

In the following contexts, we continue to describe (a) the solutions of the transition fields within stack layers via equation (70), especially for these uniform layers with isotropic materials which bring in the degenerate eigen-states, and (b) the continuum of fields conditioned at interfaces between stack layers. Consequently, a complete analysis for fields through all stacks can be fulfilled, and the associated near/far field optics can be evaluated.

#### **A.5 Eigen-system solutions**

As indicated in equation (70), the tangential fields **f**ˆ*<sup>t</sup>* within the layers proceed an eigen-system process, in which the eigen-states are independent to each other and allow individual/straightforward analyses to evaluate the transition behaviors through the layers. At the interfaces among the layers, the tangential fields **f**ˆ*<sup>t</sup>* associated with the composite phases/amplitudes of the eigen-states follow the physical continuous conditions in the laboratory framework. These characteristics lead to the necessary transform between the laboratory and eigen-system frameworks as described below. Besides, for these uniform layers with isotropic materials, especially for the incident and emitted regions, the eigen-system shows the degenerate status, and a reasonable choice of the eigen-states corresponding to the physical conditions is emphasized below. Implemented with all these, the behaviors of the tangential fields **f**ˆ*<sup>t</sup>* through all stacks layers including the in-between interfaces can be decided. The normal fields **f***n*<sup>ˆ</sup> are then obtained by equation (69), and thereby the complete light waves are understood.

#### **A.5.1 Uniform layers with isotropic materials**

For the uniform layers with isotropic materials, i.e. *ε*(*r*) is a scalar constant, the coupled-wave equation of the tangential fields **f**ˆ*<sup>t</sup>* in equation (70) can be simplified as:

$$\begin{aligned} \frac{\partial}{\partial \mathbf{\tilde{r}}} \begin{bmatrix} \vec{\tilde{e}}\_{\mathcal{X}} \\ \vec{\tilde{h}}\_{\mathcal{Y}} \\ \vec{\tilde{e}}\_{\mathcal{Y}} \\ \vec{\tilde{h}}\_{\mathcal{X}} \end{bmatrix} &= i\mathbf{C} \cdot \begin{bmatrix} \vec{\tilde{e}}\_{\mathcal{X}} \\ \vec{\tilde{h}}\_{\mathcal{Y}} \\ \vec{\tilde{e}}\_{\mathcal{Y}} \\ \vec{\tilde{h}}\_{\mathcal{X}} \end{bmatrix} \\ &= i \begin{bmatrix} 0 & \vec{n}\_{\mathcal{X}} \vec{\varepsilon}^{-1} \vec{n}\_{\mathcal{X}} - 1 & 0 & -\vec{n}\_{\mathcal{X}} \vec{\varepsilon}^{-1} \vec{n}\_{\mathcal{Y}} \\ -\vec{\varepsilon} + \vec{n}\_{\mathcal{Y}} \vec{n}\_{\mathcal{Y}} & 0 & -\vec{n}\_{\mathcal{Y}} \vec{n}\_{\mathcal{X}} & 0 \\ 0 & \vec{n}\_{\mathcal{Y}} \vec{\varepsilon}^{-1} \vec{n}\_{\mathcal{X}} & 0 & -\vec{n}\_{\mathcal{Y}} \vec{\varepsilon}^{-1} \vec{n}\_{\mathcal{Y}} + 1 \\ \vec{n}\_{\mathcal{X}} \vec{\eta}\_{\mathcal{Y}} & 0 & \vec{\varepsilon} - \vec{n}\_{\mathcal{X}} \vec{n}\_{\mathcal{X}} & 0 \end{bmatrix} \cdot \begin{bmatrix} \vec{\mathcal{E}}\_{\mathcal{X}} \\ \vec{\tilde{h}}\_{\mathcal{Y}} \\ \vec{\tilde{\varepsilon}}\_{\mathcal{Y}} \\ \vec{\tilde{\varepsilon}}\_{\mathcal{Y}} \end{bmatrix} \end{aligned} \tag{71}$$

Here, all the submatrices in **C** are diagonal and thereby the component states are independent. By straightforward calculation, its eigen-values as well as the corresponding eigen-vectors for (*g*, *h*)-order component can be obtained:

$$\begin{aligned} \textit{eigenl} & \equiv \kappa\_{\mathcal{S}^h} \\ &= \begin{bmatrix} -\mathfrak{F}\_{\mathcal{S}^h} & 0 & 0 & 0 \\ 0 & -\mathfrak{F}\_{\mathcal{S}^h} & 0 & 0 \\ 0 & 0 & \mathfrak{F}\_{\mathcal{S}^h} & 0 \\ 0 & 0 & 0 & \mathfrak{F}\_{\mathcal{S}^h} \end{bmatrix} \end{aligned} \tag{72}$$

16 Will-be-set-by-IN-TECH

In the following contexts, we continue to describe (a) the solutions of the transition fields within stack layers via equation (70), especially for these uniform layers with isotropic materials which bring in the degenerate eigen-states, and (b) the continuum of fields conditioned at interfaces between stack layers. Consequently, a complete analysis for fields through all stacks can be fulfilled, and the associated near/far field optics can be evaluated.

As indicated in equation (70), the tangential fields **f**ˆ*<sup>t</sup>* within the layers proceed an eigen-system process, in which the eigen-states are independent to each other and allow individual/straightforward analyses to evaluate the transition behaviors through the layers. At the interfaces among the layers, the tangential fields **f**ˆ*<sup>t</sup>* associated with the composite phases/amplitudes of the eigen-states follow the physical continuous conditions in the laboratory framework. These characteristics lead to the necessary transform between the laboratory and eigen-system frameworks as described below. Besides, for these uniform layers with isotropic materials, especially for the incident and emitted regions, the eigen-system shows the degenerate status, and a reasonable choice of the eigen-states corresponding to the physical conditions is emphasized below. Implemented with all these, the behaviors of the tangential fields **f**ˆ*<sup>t</sup>* through all stacks layers including the in-between interfaces can be decided. The normal fields **f***n*<sup>ˆ</sup> are then obtained by equation (69), and thereby

For the uniform layers with isotropic materials, i.e. *ε*(*r*) is a scalar constant, the coupled-wave

<sup>0</sup> *<sup>n</sup>*˜ *<sup>x</sup>ε*˜−1*n*˜ *<sup>x</sup>* <sup>−</sup> 1 0 <sup>−</sup>*n*˜ *<sup>x</sup>ε*˜

<sup>−</sup>1*n*˜ *<sup>x</sup>* <sup>0</sup> <sup>−</sup>*n*˜ *<sup>y</sup>ε*˜

−*ε*˜ + *n*˜ *yn*˜ *<sup>y</sup>* 0 −*n*˜ *yn*˜ *<sup>x</sup>* 0

*n*˜ *xn*˜ *<sup>y</sup>* 0 *ε*˜ − *n*˜ *xn*˜ *<sup>x</sup>* 0

−*ξgh* 0 00 0 −*ξgh* 0 0 0 0 *ξgh* 0 0 00 *ξgh* ⎤ ⎥ ⎥ ⎥ ⎦

Here, all the submatrices in **C** are diagonal and thereby the component states are independent. By straightforward calculation, its eigen-values as well as the corresponding eigen-vectors for

<sup>−</sup>1*n*˜ *<sup>y</sup>*

⎤ ⎥ ⎥ ⎥ ⎦ · ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(71)

(72)

<sup>−</sup>1*n*˜ *<sup>y</sup>* + 1

equation of the tangential fields **f**ˆ*<sup>t</sup>* in equation (70) can be simplified as:

0 *n*˜ *yε*˜

*eigval* ≡ *κgh*

=

⎡ ⎢ ⎢ ⎢ ⎣

**A.5 Eigen-system solutions**

the complete light waves are understood.

*∂ ∂z* ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ <sup>⎦</sup> <sup>=</sup> *<sup>i</sup>***<sup>C</sup>** ·

**A.5.1 Uniform layers with isotropic materials**

= *i*

(*g*, *h*)-order component can be obtained:

⎡ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

*ex hy ey hx* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ with *<sup>ξ</sup>gh* <sup>=</sup> �*<sup>ε</sup>* <sup>−</sup> *nyhnyh* <sup>−</sup> *nxgnxg* while the corresponding eigen-vector matrix are:

$$
\begin{split}
\begin{bmatrix}
\begin{array}{c}
\mathbf{v}'\_{g\boldsymbol{h}1} & \mathbf{v}'\_{g\boldsymbol{h}2} & \mathbf{v}'\_{g\boldsymbol{h}3} & \mathbf{v}'\_{g\boldsymbol{h}4} \\
\hline
n\_{xy}n\_{y\boldsymbol{h}} & \frac{n\_{xy}n\_{xy} - \varepsilon\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{xy}n\_{y\boldsymbol{h}}} & \frac{\mathfrak{E}\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{\boldsymbol{g}\boldsymbol{h}}n\_{y\boldsymbol{h}}} & \frac{n\_{xy}n\_{xy} - \varepsilon\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{\boldsymbol{g}\boldsymbol{h}}n\_{y\boldsymbol{h}}} \\
\hline
n\_{xy}n\_{y\boldsymbol{h}} & -\varepsilon\_{\boldsymbol{g}\boldsymbol{h}}\mathfrak{E}\_{\boldsymbol{g}\boldsymbol{h}} & \frac{n\_{y\boldsymbol{h}}\mathfrak{E}\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{\boldsymbol{g}\boldsymbol{h}}n\_{y\boldsymbol{h}}} & \frac{n\_{y\boldsymbol{h}}\mathfrak{E}\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{\boldsymbol{g}\boldsymbol{h}}n\_{y\boldsymbol{h}}} & \frac{\varepsilon\_{\boldsymbol{g}\boldsymbol{h}}\mathfrak{E}\_{\boldsymbol{g}\boldsymbol{h}}}{n\_{\boldsymbol{g}\boldsymbol{h}}n\_{y\boldsymbol{h}}} \\
0 & 1 & 0 & 1 \\
1 & 0 & 1 & 0
\end{bmatrix}
\end{bmatrix}
\tag{73}$$

Due to the degeneracy in (*κgh*,1, *κgh*,2) and (*κgh*,3, *κgh*,4), the eigenvector (**v**� *gh*1, **v**� *gh*2) as well as (**v**� *gh*3, **v**� *gh*4) can be remixed by arbitrary linear combinations. Choosing

$$m\_{\mathcal{S}^h} = \sqrt{n\_{yh}n\_{yh} + n\_{\mathcal{X}\mathcal{S}}n\_{\mathcal{X}\mathcal{S}}} \tag{74}$$

$$\mathbf{v}\_{\mathcal{g}\hbar1} = \left(\mathfrak{n}\_{\mathcal{xg}}\mathfrak{z}\_{\mathcal{g}\hbar}\mathfrak{v}'\_{\mathcal{g}\hbar1} - \mathfrak{n}\_{\mathcal{xg}}\mathfrak{v}'\_{\mathcal{g}\hbar2}\right) / m\_{\mathcal{g}\hbar} \tag{75}$$

$$\mathbf{v}\_{gh2} = \left( -\frac{\varepsilon\_{gh}\mathbf{n}\_{gh}}{\mathfrak{T}\_{gh}} \mathbf{v}'\_{gh1} + n\_{gh} \mathbf{v}'\_{gh2} \right) / m\_{gh} \tag{76}$$

$$\mathbf{v}\_{\mathcal{g}\hbar3} = \left( -n\_{\mathbf{x}\mathcal{G}} \mathfrak{z}\_{\mathcal{g}\hbar} \mathbf{v}'\_{\mathcal{g}\hbar3} - n\_{\mathbf{x}\mathcal{g}} \mathbf{v}'\_{\mathcal{g}\hbar4} \right) / m\_{\mathcal{g}\hbar} \tag{77}$$

$$\mathbf{v}\_{gh4} = \left(\frac{\varepsilon\_{gh}\mathbf{u}\_{gh}}{\mathfrak{T}\_{\mathcal{S}^{\text{fl}}}}\mathbf{v}'\_{gh3} + n\_{gh}\mathbf{v}'\_{gh4}\right) / m\_{\mathcal{g}^{\text{fl}}}\tag{78}$$

the equation (73) is then shown as:

$$
\begin{split}
\mathtt{eigenc} &= \mathtt{T}\_{\mathcal{g}h} = \begin{bmatrix}
\mathtt{v}\_{\mathcal{g}h1} & \mathtt{v}\_{\mathcal{g}h2} \ \mathtt{v}\_{\mathcal{g}h3} \ \mathtt{v}\_{\mathcal{g}h4}
\end{bmatrix} \\ &= \begin{bmatrix}
\frac{n\_{\mathcal{g}h}}{m\_{\mathcal{g}h}} & \frac{n\_{\mathcal{g}\mathcal{g}}}{m\_{\mathcal{g}h}} & \frac{n\_{\mathcal{g}h}}{m\_{\mathcal{g}h}} & \frac{n\_{\mathcal{g}\mathcal{g}}}{m\_{\mathcal{g}h}} \\
\frac{n\_{\mathcal{g}h}\mathsf{v}\_{\mathcal{g}h}}{m\_{\mathcal{g}h}} & \frac{\varepsilon\_{\mathcal{g}h}n\_{\mathcal{g}\mathcal{k}}\xi\_{\mathcal{g}h}^{-1}}{m\_{\mathcal{g}h}} & \frac{-n\_{\mathcal{g}h}\mathsf{t}\_{\mathcal{g}\mathcal{k}}}{m\_{\mathcal{g}h}} & \frac{-\varepsilon\_{\mathcal{g}h}n\_{\mathcal{k}\mathcal{k}}\xi\_{\mathcal{g}h}^{-1}}{m\_{\mathcal{g}h}} \\
\frac{-n\_{\mathcal{g}\mathcal{k}}}{m\_{\mathcal{g}h}} & \frac{n\_{\mathcal{g}h}}{m\_{\mathcal{g}h}} & \frac{-n\_{\mathcal{g}\mathcal{k}}\mathsf{t}\_{\mathcal{g}h}}{m\_{\mathcal{g}h}} & \frac{\varepsilon\_{\mathcal{g}h}n\_{\mathcal{g}\mathcal{k}}\xi\_{\mathcal{g}h}^{-1}}{m\_{\mathcal{g}h}}
\end{bmatrix} \\ \end{split} \tag{79}$$

Hence, **v***gh*<sup>1</sup> and **v***gh*<sup>2</sup> correspond to the (*g*, *h*)-order forward TE and TM (transverse electric and transverse magnetic) representations (with respect to the plane of the diffraction wave), respectively. **v***gh*<sup>3</sup> and **v***gh*<sup>4</sup> then correspond to the backward TE and TM ones. For example, with equation (69) and (79), **v***gh*<sup>1</sup> denotes the component fields:

$$\mathbf{e}\_{\mathcal{g}h} = \frac{n\_{\mathcal{y}h}}{m\_{\mathcal{g}h}} \mathbf{f} - \frac{n\_{\mathcal{xg}}}{m\_{\mathcal{g}h}} \mathbf{f} \tag{80}$$

$$\mathbf{h}\_{\mathcal{g}h} = \frac{n\_{\mathcal{xg}}\mathfrak{T}\_{\mathcal{g}h}}{m\_{\mathcal{g}h}}\mathbf{\hat{t}} + \frac{n\_{\mathcal{y}h}\mathfrak{T}\_{\mathcal{g}h}}{m\_{\mathcal{g}h}}\mathbf{\hat{t}} - m\_{\mathcal{g}h}\hat{\mathbf{k}}\tag{81}$$

**A.6 Boundary conditions**

**q**ˆ*t*,*N*+<sup>1</sup> =

� *E*− ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

or alternatively:

*<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>−</sup>

*q*,*N*+1 �*T*

 *E*<sup>+</sup> *q*,*N*+1 *M* <sup>+</sup> *q*,*N*+1 *E*− *q*,*N*+1 *M* <sup>−</sup> *q*,*N*+1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

= **T**−<sup>1</sup> *N*+1**T***<sup>a</sup>*

**q**ˆ*t*,0 ≡

� **q**<sup>+</sup> ˆ*t*,0 **q**− ˆ*t*,0 � = � **W**<sup>1</sup> **W**<sup>2</sup> **W**<sup>3</sup> **W**<sup>4</sup>

*zi*, the constriction equations can be shown as

in Three-Dimensional Periodic Liquid-Crystal Microstructures

diffraction/reflection waves can be written as:

**<sup>q</sup>**ˆ*t*,*N*+1(*zn*) = **<sup>T</sup>**−<sup>1</sup>

**T**(*a*) *<sup>i</sup>* **q**ˆ*t*,*<sup>i</sup>*

*N*+1**T**(*a*)

*N*+1**T**(*a*)

*N*+1**T**(*a*)

×(**T**(*a*)

incident region 0 and in the emitted region *N* + 1 can be obtained as:

*<sup>N</sup>* exp [*iκ<sup>a</sup>*

×... ×(**T**(*a*)

= **T**−<sup>1</sup>

= **T**−<sup>1</sup>

*<sup>N</sup>* **q**ˆ*t*,*N*(*zN*)

*iκ* (*a*)

*iκ* (*a*)

*<sup>N</sup>*−<sup>1</sup> exp �

where the first boundary is indexed as 0. Consequently, the relation between fields in the

*<sup>N</sup>* (*zN* <sup>−</sup> *zN*−1)] ...(**T***<sup>a</sup>*

Consider that the reflective field in the emitted region is zero, i.e. **q**−

<sup>ˆ</sup>*t*,*N*+<sup>1</sup> <sup>=</sup> **<sup>W</sup>**−<sup>1</sup>

**q**<sup>+</sup>

� � **q**<sup>+</sup> ˆ*t*,*N*+1 **q**− ˆ*t*,*N*+1

*<sup>N</sup>* (*zN* − *zN*−1)

*<sup>N</sup>* (*zN* − *zN*−1)

*iκ* (*a*)

> 1) <sup>−</sup><sup>1</sup> **T**<sup>0</sup>

= 0. The transmittance field in the emitted region can be obtained as:

<sup>1</sup> **<sup>q</sup>**<sup>+</sup>

�

�

*<sup>N</sup>*−<sup>1</sup> (*zN*−<sup>1</sup> <sup>−</sup> *zN*−2)

<sup>1</sup> )−1**T**0**q**ˆ*t*,0(*z*0) (88)

⎡ ⎢ ⎢ ⎢ ⎣

�

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0 ⎤ ⎥ ⎥ ⎥ <sup>⎦</sup> <sup>≡</sup> **<sup>W</sup>**−<sup>1</sup>

**<sup>q</sup>**ˆ*t*,*N*(*zN*−1)

�

⎡ ⎢ ⎢ ⎢ ⎣

= **Wq**ˆ*t*,*N*+<sup>1</sup> (90)

<sup>ˆ</sup>*t*,0 (91)

 *E*<sup>+</sup> *q*,0 *M* <sup>+</sup> *q*,0 *E*− *q*,0 *M* <sup>−</sup> *q*,0 ⎤ ⎥ ⎥ ⎥

<sup>⎦</sup> <sup>=</sup> **<sup>W</sup>**−1**q**ˆ*t*,0

<sup>ˆ</sup>*t*,*N*+<sup>1</sup> =

(89)

*<sup>N</sup>* exp �

*<sup>N</sup>* exp �

*<sup>N</sup>* )−1**T**(*a*)

Now for each layer, we have been able to independently solve the transition of electromagnetic fields in the individual layers, but the continuum of fields on interfaces has still not been included. Considering the tangential components in **f**ˆ*<sup>t</sup>* are continuous across *ith* interface at

<sup>81</sup> Electromagnetic Formalisms for Optical Propagation

(*zi*) = **<sup>T</sup>**(*a*)

Grouping this condition into the fields **q**ˆ*<sup>t</sup>* in equation (86), and introducing two virtual layers to consider the Fresnel refraction and reflection at surfaces of the media as described in the texts, a general expression for *N*−multilayered periodic structures can be obtained as in equation (13). This argument ignores the effects of multiple reflections as applied by (extended) Jones method, and similarly supplies as a easy-manipulated method. Further, an alternative process to consider the multiple reflections is described as below for references. Similarly, group the equation (87) with (86), the consecutive matrix equation with undecided

*<sup>i</sup>*+1**q**ˆ*t*,*i*+1(*zi*) (87)

along the direction **ngh** = *nxgı*<sup>ˆ</sup> + *nyh* <sup>ˆ</sup>*<sup>j</sup>* + *<sup>ξ</sup>gh* <sup>ˆ</sup> *k*. It can be seen that the characteristic fields in equations (80)-(81) associated with the eigen-solution ∝ *exp*(−*iξghz*) and constitutes the forwards TE wave. It is noted that the field amplitudes are normalized to |**e***gh*| = 1, <sup>|</sup>**h***gh*<sup>|</sup> <sup>=</sup> <sup>√</sup>*ε*, and **<sup>e</sup>***gh* · **<sup>n</sup>***gh* <sup>=</sup> **<sup>h</sup>***gh* · **<sup>n</sup>***gh* <sup>=</sup> **<sup>e</sup>***gh* · **<sup>h</sup>***gh* <sup>=</sup> 0 - that is, **<sup>n</sup>***gh*, **<sup>e</sup>***gh*, and **<sup>h</sup>***gh* are mutually perpendicular. Similarly, the remaining eigen-vectors can characterize the forwards and backwards TE/TM waves and are omitted here. In this way, a unit-amplitude incident wave then can be given as *E*<sup>+</sup> *<sup>q</sup>* = [0...1...0] *t* , *M* <sup>+</sup> *<sup>q</sup>* = 0 for TE wave, and *M* <sup>+</sup> *<sup>q</sup>* = [0...1...0] *t* , *E*<sup>+</sup> *<sup>q</sup>* = 0 for TM wave as defined below.

Considering the full components *gmin* ≤ *g* ≤ *gmax* and *hmin* ≤ *h* ≤ *hmax*, the coupled-wave equation (71) can be straightforwardly written as:

$$\begin{aligned} \frac{\partial}{\partial \boldsymbol{\Xi}} \mathbf{f}\_{\boldsymbol{f}} &= \boldsymbol{i} \mathbf{C} f\_{\boldsymbol{f}} \\ &\Rightarrow \frac{\partial}{\partial \boldsymbol{\Xi}} \mathbf{T}^{-1} \mathbf{f}\_{\boldsymbol{f}} = \boldsymbol{i} \mathbf{T}^{-1} \mathbf{C} \mathbf{T} \mathbf{T}^{-1} f\_{\boldsymbol{f}} \\ &\Rightarrow \frac{\partial}{\partial \boldsymbol{\Xi}} \mathbf{q}\_{\boldsymbol{f}} = \boldsymbol{i} \kappa \mathbf{q}\_{\boldsymbol{f}} \qquad \text{with} \quad \mathbf{f}\_{\boldsymbol{f}} = \mathbf{T} \mathbf{q}\_{\boldsymbol{f}} \end{aligned} \tag{82}$$

where

$$\mathbf{T} = \begin{bmatrix} \dot{\mathbf{n}}\_{\mathcal{Y}} & \dot{\mathbf{n}}\_{\mathcal{X}} & \dot{\mathbf{n}}\_{\mathcal{Y}} & \dot{\mathbf{n}}\_{\mathcal{X}} \\ \dot{\mathbf{n}}\_{\mathcal{Y}}\tilde{\boldsymbol{\xi}} & \varepsilon \dot{\mathbf{n}}\_{\mathcal{X}}\tilde{\boldsymbol{\xi}}^{-1} & -\dot{\mathbf{n}}\_{\mathcal{Y}}\tilde{\boldsymbol{\xi}} - \varepsilon \dot{\mathbf{n}}\_{\mathcal{X}}\tilde{\boldsymbol{\xi}}^{-1} \\ -\dot{\mathbf{n}}\_{\mathcal{X}} & \dot{\mathbf{n}}\_{\mathcal{Y}} & -\dot{\mathbf{n}}\_{\mathcal{X}} & \dot{\mathbf{n}}\_{\mathcal{Y}} \\ \dot{\mathbf{n}}\_{\mathcal{X}}\tilde{\boldsymbol{\xi}} & -\varepsilon \dot{\mathbf{n}}\_{\mathcal{Y}}\tilde{\boldsymbol{\xi}}^{-1} - \dot{\mathbf{n}}\_{\mathcal{X}}\tilde{\boldsymbol{\xi}} & \varepsilon \dot{\mathbf{n}}\_{\mathcal{Y}}\tilde{\boldsymbol{\xi}}^{-1} \end{bmatrix}, \quad \mathbf{q}\_{\mathcal{I}} = \begin{bmatrix} \vec{E}\_{q}^{+} \\ \vec{M}\_{q}^{+} \\ \vec{E}\_{q}^{-} \\ \vec{M}\_{q}^{-} \end{bmatrix} \tag{83}$$

Note that ˙**n***<sup>y</sup>* and ˙**n***<sup>x</sup>* are the *NgNh* <sup>×</sup> *NgNh* diagonal matrices with diagonal elements *nyh mgh* and *nxg mgh* respectively. *<sup>ξ</sup>*−<sup>1</sup> is the matrix with elements 1/*ξgh*, not the inverse of the matrix *<sup>ξ</sup>*. Moreover, *E*<sup>+</sup> *<sup>q</sup>* and *M* <sup>+</sup> *<sup>q</sup>* ( *E*− *<sup>q</sup>* and *M* <sup>−</sup> *<sup>q</sup>* ) correspond to the physical forward (backward) *TE* and *TM* waves, respectively. The transition of fields **q**ˆ*<sup>t</sup>* within the considered layer are now solved as:

$$\mathbf{q}\_{\rm f}(\overline{z}) = e \exp\left[i\kappa \left(\overline{z} - \overline{z}\_{0}\right)\right] \mathbf{q}\_{\rm f}(\overline{z}\_{0}) \tag{84}$$

#### **A.5.2 Periodic-structured layers with isotropic/birefringent materials**

For the in-between periodic layers, the transition equations of tangential fields **f**ˆ*<sup>t</sup>* in equation (70) can be generally expressed as:

$$\frac{\partial}{\partial \mathbf{2}} \mathbf{q}\_{\mathbf{f}} = i \kappa^{(a)} \mathbf{q}\_{\mathbf{f}} \qquad \text{with} \quad \mathbf{f}\_{\mathbf{f}} = \mathbf{T}^{(a)} \mathbf{q}\_{\mathbf{f}} \tag{85}$$

with the transition of **q**ˆ*<sup>t</sup>*

$$\mathbf{q}\_f(\overline{\mathbf{z}}) = \exp\left[i\kappa^{(a)}\left(\overline{\mathbf{z}} - \overline{\mathbf{z}}\_0\right)\right] \mathbf{q}\_f(\overline{\mathbf{z}}\_0) \tag{86}$$

Here, **T**(*a*) is the eigen-vector matrix of **G** of equation (70) with column eigen-vectors, and *κ*(*a*) is the corresponding diagonal eigen-value matrix.

#### **A.6 Boundary conditions**

18 Will-be-set-by-IN-TECH

in equations (80)-(81) associated with the eigen-solution ∝ *exp*(−*iξghz*) and constitutes the forwards TE wave. It is noted that the field amplitudes are normalized to |**e***gh*| = 1, <sup>|</sup>**h***gh*<sup>|</sup> <sup>=</sup> <sup>√</sup>*ε*, and **<sup>e</sup>***gh* · **<sup>n</sup>***gh* <sup>=</sup> **<sup>h</sup>***gh* · **<sup>n</sup>***gh* <sup>=</sup> **<sup>e</sup>***gh* · **<sup>h</sup>***gh* <sup>=</sup> 0 - that is, **<sup>n</sup>***gh*, **<sup>e</sup>***gh*, and **<sup>h</sup>***gh* are mutually perpendicular. Similarly, the remaining eigen-vectors can characterize the forwards and backwards TE/TM waves and are omitted here. In this way, a unit-amplitude incident

Considering the full components *gmin* ≤ *g* ≤ *gmax* and *hmin* ≤ *h* ≤ *hmax*, the coupled-wave

**<sup>T</sup>**−1**f**ˆ*<sup>t</sup>* <sup>=</sup> *<sup>i</sup>***T**−1**CTT**−<sup>1</sup> *<sup>f</sup>*ˆ*<sup>t</sup>*

*k*. It can be seen that the characteristic fields

*<sup>q</sup>* = [0...1...0]

*t* , *E*<sup>+</sup> *<sup>q</sup>* = 0

(83)

*mgh*

*<sup>q</sup>* = 0 for TE wave, and *M* <sup>+</sup>

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

, **q**ˆ*<sup>t</sup>* =

*<sup>q</sup>* ) correspond to the physical forward (backward) *TE* and

**q**ˆ*t*(*z*) = *exp* [*iκ* (*z* − *z*0)] **q**ˆ*t*(*z*0) (84)

**q**ˆ*<sup>t</sup>* = *iκ***q**ˆ*<sup>t</sup> with* **f**ˆ*<sup>t</sup>* = **Tq**ˆ*<sup>t</sup>* (82)

⎡ ⎢ ⎢ ⎢ ⎣

 *E*<sup>+</sup> *q M* <sup>+</sup> *q E*− *q M* <sup>−</sup> *q*

⎤ ⎥ ⎥ ⎥ ⎦

**q**ˆ*<sup>t</sup>* (85)

**q**ˆ*t*(*z*0) (86)

along the direction **ngh** = *nxgı*<sup>ˆ</sup> + *nyh* <sup>ˆ</sup>*<sup>j</sup>* + *<sup>ξ</sup>gh* <sup>ˆ</sup>

*E*<sup>+</sup>

*∂ ∂z*

equation (71) can be straightforwardly written as:

**T** =

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *<sup>q</sup>* = [0...1...0]

**f**ˆ*<sup>t</sup>* = *i***C** *f*ˆ*<sup>t</sup>*

⇒ *∂ ∂z*

⇒ *∂ ∂z*

**n**˙ *y* **n**˙ *x* **n**˙ *y* **n**˙ *x* **<sup>n</sup>**˙ *<sup>y</sup>ξ ε***n**˙ *<sup>x</sup>ξ*−<sup>1</sup> <sup>−</sup>**n**˙ *<sup>y</sup><sup>ξ</sup>* <sup>−</sup>*ε***n**˙ *<sup>x</sup>ξ*−<sup>1</sup> −**n**˙ *<sup>x</sup>* **n**˙ *<sup>y</sup>* −**n**˙ *<sup>x</sup>* **n**˙ *<sup>y</sup>* **<sup>n</sup>**˙ *<sup>x</sup><sup>ξ</sup>* <sup>−</sup>*ε***n**˙ *<sup>y</sup>ξ*−<sup>1</sup> <sup>−</sup>**n**˙ *<sup>x</sup>ξ ε***n**˙ *<sup>y</sup>ξ*−<sup>1</sup>

Note that ˙**n***<sup>y</sup>* and ˙**n***<sup>x</sup>* are the *NgNh* <sup>×</sup> *NgNh* diagonal matrices with diagonal elements *nyh*

*TM* waves, respectively. The transition of fields **q**ˆ*<sup>t</sup>* within the considered layer are now solved

For the in-between periodic layers, the transition equations of tangential fields **f**ˆ*<sup>t</sup>* in equation

**<sup>q</sup>**ˆ*<sup>t</sup> with* **<sup>f</sup>**ˆ*<sup>t</sup>* <sup>=</sup> **<sup>T</sup>**(*a*)

�

*<sup>i</sup>κ*(*a*) (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*0)

Here, **T**(*a*) is the eigen-vector matrix of **G** of equation (70) with column eigen-vectors, and

*mgh* respectively. *<sup>ξ</sup>*−<sup>1</sup> is the matrix with elements 1/*ξgh*, not the inverse of the matrix *<sup>ξ</sup>*.

*t* , *M* <sup>+</sup>

wave then can be given as

where

and *nxg*

as:

Moreover,

*E*<sup>+</sup>

with the transition of **q**ˆ*<sup>t</sup>*

*<sup>q</sup>* and *M* <sup>+</sup>

(70) can be generally expressed as:

*<sup>q</sup>* ( *E*−

*<sup>q</sup>* and *M* <sup>−</sup>

*∂ ∂z*

*κ*(*a*) is the corresponding diagonal eigen-value matrix.

**A.5.2 Periodic-structured layers with isotropic/birefringent materials**

**<sup>q</sup>**ˆ*<sup>t</sup>* <sup>=</sup> *<sup>i</sup>κ*(*a*)

**q**ˆ*t*(*z*) = *exp*

�

for TM wave as defined below.

Now for each layer, we have been able to independently solve the transition of electromagnetic fields in the individual layers, but the continuum of fields on interfaces has still not been included. Considering the tangential components in **f**ˆ*<sup>t</sup>* are continuous across *ith* interface at *zi*, the constriction equations can be shown as

$$\mathbf{T}\_{i}^{(a)}\mathbf{q}\_{\hat{t},i}(\overline{z}\_{i}) = \mathbf{T}\_{i+1}^{(a)}\mathbf{q}\_{\hat{t},i+1}(\overline{z}\_{i})\tag{87}$$

Grouping this condition into the fields **q**ˆ*<sup>t</sup>* in equation (86), and introducing two virtual layers to consider the Fresnel refraction and reflection at surfaces of the media as described in the texts, a general expression for *N*−multilayered periodic structures can be obtained as in equation (13). This argument ignores the effects of multiple reflections as applied by (extended) Jones method, and similarly supplies as a easy-manipulated method. Further, an alternative process to consider the multiple reflections is described as below for references.

Similarly, group the equation (87) with (86), the consecutive matrix equation with undecided diffraction/reflection waves can be written as:

$$\begin{split} \mathbf{q}\_{\hat{\mathbf{f}},N+1}(\boldsymbol{\Xi}\_{n}) &= \mathbf{T}\_{N+1}^{-1} \mathbf{T}\_{N}^{(a)} \mathbf{q}\_{\hat{\mathbf{f}},N}(\boldsymbol{\Xi}\_{N}) \\ &= \mathbf{T}\_{N+1}^{-1} \mathbf{T}\_{N}^{(a)} \exp\left[i\kappa\_{N}^{(a)} (\boldsymbol{\Xi}\_{N} - \boldsymbol{\Xi}\_{N-1})\right] \mathbf{q}\_{\hat{\mathbf{f}},N}(\boldsymbol{\Xi}\_{N-1}) \\ &= \mathbf{T}\_{N+1}^{-1} \mathbf{T}\_{N}^{(a)} \exp\left[i\kappa\_{N}^{(a)} (\boldsymbol{\Xi}\_{N} - \boldsymbol{\Xi}\_{N-1})\right] \\ &\quad \times (\mathbf{T}\_{N}^{(a)})^{-1} \mathbf{T}\_{N-1}^{(a)} \exp\left[i\kappa\_{N-1}^{(a)} (\boldsymbol{\Xi}\_{N-1} - \boldsymbol{\Xi}\_{N-2})\right] \\ &\quad \times \dots \\ &\quad \times (\mathbf{T}\_{1}^{(a)})^{-1} \mathbf{T}\_{0} \mathbf{q}\_{\hat{\mathbf{f}},0}(\boldsymbol{\Xi}\_{0}) \end{split} \tag{88}$$

where the first boundary is indexed as 0. Consequently, the relation between fields in the incident region 0 and in the emitted region *N* + 1 can be obtained as:

$$\mathbf{q}\_{l,N+1} = \begin{bmatrix} \vec{E}\_{g,N+1}^{+} \\ \vec{M}\_{g,N+1}^{+} \\ \vec{E}\_{g,N+1}^{-} \\ \vec{M}\_{g,N+1}^{-} \end{bmatrix} = \mathbf{T}\_{N+1}^{-1} \mathbf{T}\_{N}^{q} \exp\left[i \mathbf{x}\_{N}^{a} \left(\Xi\_{N} - \Xi\_{N-1}\right)\right] \dots \left(\mathbf{T}\_{1}^{q}\right)^{-1} \mathbf{T}\_{0} \begin{bmatrix} \vec{E}\_{g,0}^{+} \\ \vec{M}\_{g,0}^{+} \\ \vec{E}\_{g,0}^{-} \\ \vec{M}\_{g,0}^{-} \end{bmatrix} \equiv \mathbf{W}^{-1} \begin{bmatrix} \vec{E}\_{g,0}^{+} \\ \vec{M}\_{g,0}^{+} \\ \vec{E}\_{g,0}^{-} \\ \vec{M}\_{g,0}^{-} \end{bmatrix} = \mathbf{W}^{-1} \mathbf{q}\_{l,0} \tag{89}$$

or alternatively:

$$\mathbf{q}\_{\hat{t},0} \equiv \begin{bmatrix} \mathbf{q}\_{\hat{t},0}^{+} \\ \mathbf{q}\_{\hat{t},0}^{-} \end{bmatrix} = \begin{bmatrix} \mathbf{W}\_1 \ \mathbf{W}\_2 \\ \mathbf{W}\_3 \ \mathbf{W}\_4 \end{bmatrix} \begin{bmatrix} \mathbf{q}\_{\hat{t},N+1}^{+} \\ \mathbf{q}\_{\hat{t},N+1}^{-} \end{bmatrix} = \mathbf{W} \mathbf{q}\_{\hat{t},N+1} \tag{90}$$

Consider that the reflective field in the emitted region is zero, i.e. **q**− <sup>ˆ</sup>*t*,*N*+<sup>1</sup> = � *E*− *<sup>q</sup>*,*N*+<sup>1</sup> *<sup>M</sup>* <sup>−</sup> *q*,*N*+1 �*T* = 0. The transmittance field in the emitted region can be obtained as:

$$\mathbf{q}\_{\hat{\mathbf{f}},N+1}^{+} = \mathbf{W}\_1^{-1} \mathbf{q}\_{\hat{\mathbf{f}},0}^{+} \tag{91}$$

*ε*yzgh=InverseFourier[(ne∧2 − no∧2)\*Sin[*θ*o]Cos[*θ*o]Sin[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy];

<sup>83</sup> Electromagnetic Formalisms for Optical Propagation

*ε*xx = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}];*ε*xy = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}]; *ε*xz = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}];*ε*yy = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}]; *ε*yz = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}];*ε*zz = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}];

*ε*xx[[*i*, *j*]] = *ε*xxgh[[gp, hp]];*ε*xy[[*i*, *j*]] = *ε*xygh[[gp, hp]];*ε*xz[[*i*, *j*]] = *ε*xzgh[[gp, hp]]; *ε*yy[[*i*, *j*]] = *ε*yygh[[gp, hp]];*ε*yz[[*i*, *j*]] = *ε*yzgh[[gp, hp]];*ε*zz[[*i*, *j*]] = *ε*zzgh[[gp, hp]];

G11 = Dot[nx,*ε*zzinv,*ε*xz]; G12 = Dot[nx,*ε*zzinv, nx] − IdentityMatrix[Ng\*Nh];

G33 = Dot[ny,*ε*zzinv,*ε*yz]; G34 = −Dot[ny,*ε*zzinv, ny] + IdentityMatrix[Ng\*Nh]; G41 = −Dot[*ε*yz,*ε*zzinv,*ε*xz] + *ε*xy + Dot[nx, ny]; G42 = −Dot[*ε*yz,*ε*zzinv, nx]; G43 = −Dot[*ε*yz,*ε*zzinv,*ε*yz] + *ε*yy − Dot[nx, nx]; G44 = Dot[*ε*yz,*ε*zzinv, ny];

*κ*a = Dot[Tainv, *G*, Ta]; (\*eigen − value matrix corresponding to the arrangement of Ta\*)

(\*Calculate the matrixes related to incidnet and emitted air regions, *i*.*e*. nI = nE = 1\*)

G21 = Dot[*ε*xz, *ε*zzinv,*ε*xz] − *ε*xx + Dot[ny, ny]; G22 = Dot[*ε*xz,*ε*zzinv, nx]; G23 = Dot[*ε*xz, *ε*zzinv,*ε*yz] − *ε*xy − Dot[ny, nx]; G24 = −Dot[*ε*xz,*ε*zzinv, ny];

G1i = Join[G11, G12, G13, G14, 2]; G2i = Join[G21, G22, G23, G24, 2]; G3i = Join[G31, G32, G33, G34, 2]; G4i = Join[G41, G42, G43, G44, 2];

Ta = Transpose[Eigenvectors[*G*]]; (\*eigen − vecotr matrix\*) Tainv = Inverse[Ta]; (\*inverse of the eigen − vecotr matrix \*)

nxd = DiagonalMatrix[Table[nx[[*i*, *i*]]/*m*[[*i*,*i*]], {*i*, Ng\*Nh}]]; nyd = DiagonalMatrix[Table[ny[[*i*, *i*]]/*m*[[*i*,*i*]], {*i*, Ng\*Nh}]];

T21 = Dot[nyd, *ξ*]; T22 = nI∧2\*Dot[nxd, *ξ*inv]; T23 = −Dot[nyd, *ξ*];

T41 = Dot[nxd, *ξ*]; T42 = −nI∧2\*Dot[nyd, *ξ*inv]; T43 = −Dot[nxd, *ξ*];

T1i = Join[T11, T12, T13, T14, 2]; T2i = Join[T21, T22, T23, T24, 2]; T3i = Join[T31, T32, T33, T34, 2]; T4i = Join[T41, T42, T43, T44, 2];

T11 = nyd; T12 = nxd; T13 = nyd; T14 = nxd;

T31 = −nxd; T32 = nyd; T33 = −nxd; T34 = nyd;

Ti = Join[T1i, T2i, T3i, T4i]; (\* transform matrix Ti\*)

Tiinv = Inverse[Ti]; (\* inverse of the transform matrix Ti \*)

∧2]/Sqrt[GridNx]/Sqrt[GridNy];

*ε*zzgh=InverseFourier[no∧2 + (ne∧2 − no∧2)\*Cos[*θ*o]

in Three-Dimensional Periodic Liquid-Crystal Microstructures

*g* = gindx[[*i*]] − gindx[[*j*]]; *h* = hindx[[*i*]] − hindx[[*j*]]; gp = If[*g* ≥ 0, *g* = *g* + 1, *g* = *g* + GridNx + 1]; (\* follow arrangements of components in *ε*ijgh \*) hp = If[*h* ≥ 0, *h* = *h* + 1, *h* = *h* + GridNy + 1];

(\* Define the matrix *ε*ij with element *ε*ijgh \*)

*ε*zzinv = Table[0, {*i*, Ng\*Nh}, {*j*, Ng\*Nh}]; For[*i* = 1, *i* ≤ Ng\*Nh, For[*j* = 1, *j* ≤ Ng\*Nh,

(\* Calculate matrix G for the single LC layer\*)

G13 = Dot[nx,*ε*zzinv,*ε*yz]; G14 = −Dot[nx,*ε*zzinv, ny];

G31 = Dot[ny,*ε*zzinv,*ε*xz]; G32 = Dot[ny,*ε*zzinv, nx];

*j*++; ]; *i*++; ];

*ε*zzinv = Inverse[*ε*zz];

*G* = Join[G1i, G2i, G3i, G4i];

T24 = −nI∧2\*Dot[nxd, *ξ*inv];

T44 = nI∧2\*Dot[nyd, *ξ*inv];

and the reflective field in the incident region is:

$$\mathbf{q}\_{\hat{\mathbf{f}},0}^{-} = \mathbf{W}\_{\mathbf{3}} \mathbf{W}\_{1}^{-1} \mathbf{q}\_{\hat{\mathbf{f}},0}^{+} \tag{92}$$

#### **A.7 Diffraction efficiency**

To evaluate the diffraction efficiency with the obtained **q**± <sup>ˆ</sup>*<sup>t</sup>* , the *<sup>x</sup>*, *<sup>y</sup>*, and *<sup>z</sup>* components of the transmittance/reflection fields of the diffraction order (*g*, *h*) can be calculated by equations 19 (or 82, 83) and 69 for emitted/incident regions, and thereby the standard definitions of diffraction efficiency can be followed. Note that the incident fields should be excluded when calculating the reflection fields in the incident region.

#### **B. Program codes of Wolfram Mathematica for Coupling Matrix Method**

In this appendix, the program codes of Wolfram Mathematica for the (numerical) study case in the previous section are added as follows. It could be able to do the simulations by copy and paste the codes, while few characters may need to be adjusted, e.g., the superscript of *W*� (*W*��) and the power symbol on *no*∧2 (*ne*∧2).

(\*Initialize one − period LC profiles (*θ*o, *φ*o) for single LC layer\*) dx = 0.1; dy = dx; (\*um/grid; grid interval \*) GridNx = 100; GridNy = 100; (\*grid num. in *x* and *y* \*) Λx = GridNx\*dx; Λy = GridNx\*dy; (\* unit cell \*) *θ*o = Table[*π*\**i*/GridNx, {*i*, GridNx}, {*j*, GridNy}]; *φ*o = Table[*π*/2.0, {*i*, GridNx}, {*j*, GridNy}]; dz = 2.0; (\*um; the thickness of the LC layer \*)

(\*Define optical − related parameters\*) nI = 1.0; nE = 1.0; *θ* = 0.001; *φ* = 0.0; *λ* = 0.55; ne = 1.5; no = 1.6; grng = 1; hrng = 1; (\* − grng ≤ *g* ≤ grng; −hrng ≤ *h* ≤ hrng\*) Ng = 2\*grng + 1; Nh = 2\*hrng + 1; (\*Note Ng < GridNx, Nh < GridNy\*)

(\*Initialize relevant wave − vector matrixes related to nxg, nyh, respectively\*) gindx = Table[Floor[(*i* − 1.0)/Nh] − grng, {*i*, Ng\*Nh}]; (\* g sequence in ei or hi fields \*) hindx = Table[Mod[(*i* − 1), Nh] − hrng, {*i*, Ng\*Nh}]; (\* h sequence in ei or hi fields \*) nx = DiagonalMatrix[Table[nI\*Sin[*θ*]\*Cos[*φ*] − gindx[[*i*]]\**λ*/Λx, {*i*, Ng\*Nh}]]; ny = DiagonalMatrix[Table[nI\*Sin[*θ*]\*Sin[*φ*] − hindx[[*i*]]\**λ*/Λy, {*i*, Ng\*Nh}]]; *m* = DiagonalMatrix[Table[Sqrt[nx[[*i*, *i*]]∧2 + ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]]; *ξ* = DiagonalMatrix[Table[Sqrt[nI∧2 − nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]]; *ξ*inv = DiagonalMatrix[Table[1.0/Sqrt[nI∧2 − nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]];

(\*Calculate *ε*ijgh by Fourier transform of *ε*ij(*x*, *y*; *z*) for the single LC layer\*) *ε*xxgh=InverseFourier[no∧2+(ne∧2−no∧2)\*Sin[*θ*o] ∧2\*Cos[*φ*o] ∧2]/Sqrt[GridNx]/Sqrt[GridNy]; *ε*xygh=InverseFourier[(ne∧2−no∧2)\*Sin[*θ*o] ∧2\*Sin[*φ*o]Cos[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy]; *ε*xzgh=InverseFourier[(ne∧2−no∧2)\*Sin[*θ*o]Cos[*θ*o]Cos[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy]; *ε*yygh=InverseFourier[no∧2+(ne∧2−no∧2)\*Sin[*θ*o] ∧2\*Sin[*φ*o] ∧2]/Sqrt[GridNx]/Sqrt[GridNy];

20 Will-be-set-by-IN-TECH

<sup>ˆ</sup>*t*,0 <sup>=</sup> **<sup>W</sup>**3**W**−<sup>1</sup>

transmittance/reflection fields of the diffraction order (*g*, *h*) can be calculated by equations 19 (or 82, 83) and 69 for emitted/incident regions, and thereby the standard definitions of diffraction efficiency can be followed. Note that the incident fields should be excluded when

In this appendix, the program codes of Wolfram Mathematica for the (numerical) study case in the previous section are added as follows. It could be able to do the simulations by copy and paste the codes, while few characters may need to be adjusted, e.g., the superscript of *W*�

<sup>1</sup> **<sup>q</sup>**<sup>+</sup>

<sup>ˆ</sup>*t*,0 (92)

<sup>ˆ</sup>*<sup>t</sup>* , the *<sup>x</sup>*, *<sup>y</sup>*, and *<sup>z</sup>* components of the

∧2]/Sqrt[GridNx]/Sqrt[GridNy];

∧2]/Sqrt[GridNx]/Sqrt[GridNy];

∧2\*Sin[*φ*o]Cos[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy];

**q**−

**B. Program codes of Wolfram Mathematica for Coupling Matrix Method**

and the reflective field in the incident region is:

To evaluate the diffraction efficiency with the obtained **q**±

(\*Initialize one − period LC profiles (*θ*o, *φ*o) for single LC layer\*)

grng = 1; hrng = 1; (\* − grng ≤ *g* ≤ grng; −hrng ≤ *h* ≤ hrng\*)

Ng = 2\*grng + 1; Nh = 2\*hrng + 1; (\*Note Ng < GridNx, Nh < GridNy\*)

(\*Calculate *ε*ijgh by Fourier transform of *ε*ij(*x*, *y*; *z*) for the single LC layer\*)

*ε*xxgh=InverseFourier[no∧2+(ne∧2−no∧2)\*Sin[*θ*o]

*ε*yygh=InverseFourier[no∧2+(ne∧2−no∧2)\*Sin[*θ*o]

*ε*xygh=InverseFourier[(ne∧2−no∧2)\*Sin[*θ*o]

(\*Initialize relevant wave − vector matrixes related to nxg, nyh, respectively\*)

gindx = Table[Floor[(*i* − 1.0)/Nh] − grng, {*i*, Ng\*Nh}]; (\* g sequence in ei or hi fields \*) hindx = Table[Mod[(*i* − 1), Nh] − hrng, {*i*, Ng\*Nh}]; (\* h sequence in ei or hi fields \*) nx = DiagonalMatrix[Table[nI\*Sin[*θ*]\*Cos[*φ*] − gindx[[*i*]]\**λ*/Λx, {*i*, Ng\*Nh}]]; ny = DiagonalMatrix[Table[nI\*Sin[*θ*]\*Sin[*φ*] − hindx[[*i*]]\**λ*/Λy, {*i*, Ng\*Nh}]]; *m* = DiagonalMatrix[Table[Sqrt[nx[[*i*, *i*]]∧2 + ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]]; *ξ* = DiagonalMatrix[Table[Sqrt[nI∧2 − nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]]; *ξ*inv = DiagonalMatrix[Table[1.0/Sqrt[nI∧2 − nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]];

*ε*xzgh=InverseFourier[(ne∧2−no∧2)\*Sin[*θ*o]Cos[*θ*o]Cos[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy];

∧2\*Cos[*φ*o]

∧2\*Sin[*φ*o]

calculating the reflection fields in the incident region.

(*W*��) and the power symbol on *no*∧2 (*ne*∧2).

dx = 0.1; dy = dx; (\*um/grid; grid interval \*)

(\*Define optical − related parameters\*) nI = 1.0; nE = 1.0; *θ* = 0.001; *φ* = 0.0; *λ* = 0.55;

ne = 1.5; no = 1.6;

GridNx = 100; GridNy = 100; (\*grid num. in *x* and *y* \*) Λx = GridNx\*dx; Λy = GridNx\*dy; (\* unit cell \*) *θ*o = Table[*π*\**i*/GridNx, {*i*, GridNx}, {*j*, GridNy}]; *φ*o = Table[*π*/2.0, {*i*, GridNx}, {*j*, GridNy}]; dz = 2.0; (\*um; the thickness of the LC layer \*)

**A.7 Diffraction efficiency**

*ε*yzgh=InverseFourier[(ne∧2 − no∧2)\*Sin[*θ*o]Cos[*θ*o]Sin[*φ*o]]/Sqrt[GridNx]/Sqrt[GridNy]; *ε*zzgh=InverseFourier[no∧2 + (ne∧2 − no∧2)\*Cos[*θ*o] ∧2]/Sqrt[GridNx]/Sqrt[GridNy];

```
(* Define the matrix εij with element εijgh *)
εxx = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];εxy = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];
εxz = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];εyy = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];
εyz = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];εzz = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];
εzzinv = Table[0, {i, Ng*Nh}, {j, Ng*Nh}];
For[i = 1, i ≤ Ng*Nh, For[j = 1, j ≤ Ng*Nh,
g = gindx[[i]] − gindx[[j]]; h = hindx[[i]] − hindx[[j]];
gp = If[g ≥ 0, g = g + 1, g = g + GridNx + 1];
(* follow arrangements of components in εijgh *)
hp = If[h ≥ 0, h = h + 1, h = h + GridNy + 1];
εxx[[i, j]] = εxxgh[[gp, hp]];εxy[[i, j]] = εxygh[[gp, hp]];εxz[[i, j]] = εxzgh[[gp, hp]];
εyy[[i, j]] = εyygh[[gp, hp]];εyz[[i, j]] = εyzgh[[gp, hp]];εzz[[i, j]] = εzzgh[[gp, hp]];
j++; ]; i++; ];
εzzinv = Inverse[εzz];
(* Calculate matrix G for the single LC layer*)
G11 = Dot[nx,εzzinv,εxz]; G12 = Dot[nx,εzzinv, nx] − IdentityMatrix[Ng*Nh];
G13 = Dot[nx,εzzinv,εyz]; G14 = −Dot[nx,εzzinv, ny];
G21 = Dot[εxz, εzzinv,εxz] − εxx + Dot[ny, ny]; G22 = Dot[εxz,εzzinv, nx];
G23 = Dot[εxz, εzzinv,εyz] − εxy − Dot[ny, nx]; G24 = −Dot[εxz,εzzinv, ny];
G31 = Dot[ny,εzzinv,εxz]; G32 = Dot[ny,εzzinv, nx];
G33 = Dot[ny,εzzinv,εyz]; G34 = −Dot[ny,εzzinv, ny] + IdentityMatrix[Ng*Nh];
G41 = −Dot[εyz,εzzinv,εxz] + εxy + Dot[nx, ny]; G42 = −Dot[εyz,εzzinv, nx];
G43 = −Dot[εyz,εzzinv,εyz] + εyy − Dot[nx, nx]; G44 = Dot[εyz,εzzinv, ny];
G1i = Join[G11, G12, G13, G14, 2]; G2i = Join[G21, G22, G23, G24, 2];
G3i = Join[G31, G32, G33, G34, 2]; G4i = Join[G41, G42, G43, G44, 2];
G = Join[G1i, G2i, G3i, G4i];
Ta = Transpose[Eigenvectors[G]]; (*eigen − vecotr matrix*)
Tainv = Inverse[Ta]; (*inverse of the eigen − vecotr matrix *)
κa = Dot[Tainv, G, Ta]; (*eigen − value matrix corresponding to the arrangement of Ta*)
(*Calculate the matrixes related to incidnet and emitted air regions, i.e. nI = nE = 1*)
nxd = DiagonalMatrix[Table[nx[[i, i]]/m[[i,i]], {i, Ng*Nh}]];
nyd = DiagonalMatrix[Table[ny[[i, i]]/m[[i,i]], {i, Ng*Nh}]];
T11 = nyd; T12 = nxd; T13 = nyd; T14 = nxd;
T21 = Dot[nyd, ξ]; T22 = nI∧2*Dot[nxd, ξinv]; T23 = −Dot[nyd, ξ];
T24 = −nI∧2*Dot[nxd, ξinv];
T31 = −nxd; T32 = nyd; T33 = −nxd; T34 = nyd;
T41 = Dot[nxd, ξ]; T42 = −nI∧2*Dot[nyd, ξinv]; T43 = −Dot[nxd, ξ];
T44 = nI∧2*Dot[nyd, ξinv];
T1i = Join[T11, T12, T13, T14, 2]; T2i = Join[T21, T22, T23, T24, 2];
T3i = Join[T31, T32, T33, T34, 2]; T4i = Join[T41, T42, T43, T44, 2];
Ti = Join[T1i, T2i, T3i, T4i]; (* transform matrix Ti*)
Tiinv = Inverse[Ti]; (* inverse of the transform matrix Ti *)
```
(\* Calculate the diffraction and reflection fields \*)

in Three-Dimensional Periodic Liquid-Crystal Microstructures

Print["TE mode without multi-reflections"];

Print[gindx[[*i*]], ", ", hindx[[*i*]], ", ", Abs[diff2[[*i*]]]];*i*++; ];

= 1.0; (\* ignore backward \*) diff2 = Dot[Sext, S1, Sent, indfd2];

For[*i* = 1, *i* ≤ Ng\*Nh,

502-510.

031114-3.

1332-1341.

pp. 1201-1212.

March 2007, pp. 75-90.

**5. References**

indfd2 = Table[0, {*i*, 4\*Ng\*Nh}]; indfd2[[Round[(Ng\*Nh + 1)/2]]]

(\*Print diffraction and reflection fields as well as the corresponding *g*, *h* orders\*)

Berreman, D. W. (1972). Optics in Stratified and Anisotropic Media: 4-by-4-Matrix

<sup>85</sup> Electromagnetic Formalisms for Optical Propagation

Blinov, L. M.; Cipparrone, G.; Pagliusi, P.; Lazarev, V. V. & Palto, S. P. (2006). Mirrorless lasing

Blinov, L. M.; Lazarev, V. V.; Palto, S. P.; Cipparrone, G.; Mazzulla, A. & Pagliusi, P. (2007).

Galatola, P; Oldano, C & Kumar, P. B. S. (1994). Symmetry properties of anisotropic dielectric

Glytsis, E. N. & Gaylord, T. K. (1987), Rigorous three-dimensional coupled-wave diffraction

Ho, I. L.; Chang, Y. C.; Huang, C. H.& Li, W. Y. (2011), A detailed derivation of

Kriezisa, E. E. & Elston, S. J. (1999). A wide angle beam propagation method for the analysis of

Kriezis, E. E.; Filippov, S. K. & Elston, S. J. (2000). Light propagation in domain walls

microstructures, *Liquid Crystals*, Vol. 38, No. 2, February 2011, 241 ˛aV252. Jones, R. C. (1941). A new calculus for the treatment of optical systems, *Journal of the Optical*

*America A*, Vol. 4, Iss. 11, November 1987, pp. 2061-2080.

*Society of America*, Vol. 31, Iss. 7, July 1941, pp. 488 ˛aV493.

Issue 25, September 2002 , pp. 5346-5356.

Formulation, *Journal of the Optical Society of America*, Vol. 62, Iss. 4, April 1972, pp.

from nematic liquid crystals in the plane waveguide geometry without refractive index or gain modulation, *Applied Physics Letters*, Vol. 89, Iss. 3, July 2006, pp.

Electric field tuning a spectrum of nematic liquid crystal lasing with the use of a periodic shadow mask, *Journal of Nonlinear Optical Physics & Materials*, Vol. 16, Iss. 1,

gratings, *Journal of the Optical Society of America A*, Vol. 11, Iss. 4, April 1994, pp.

analysis of single and cascaded anisotropic gratings, *Journal of the Optical Society of*

rigorous coupled wave algorithms for three-dimensional periodic liquid-crystal

tilted nematic liquid crystal structures, *Journal of Modern Optics*, Vol. 46, Iss. 8, 1999,

in ferroelectric liquid crystal devices by the finite-difference time-domain method, *Journal of Optics A: Pure and Applied Optics*, Vol. 2, No. 1, January 2000, pp. 27-33. Kriezis, E. E. & Elston, S. J. (2000). Wide-angle beam propagation method for liquid-crystal device calculations, *Applied Optics*, Vol. 39, Iss. 31, November 2000, pp.5707-5714. Kriezis, E. E.; Newton, C. J. P.; Spiller, T. P. & Elston, S. J. (2002). Three-dimensional simulations

of light propagation in periodic liquid-crystal microstructures, *Applied Optics*, Vol. 41,

(\*Solution(1) : solve diffractions and reflections with multi − reflections\*) (\*set incident plane wave indfd,*e*.*g*. set the value of the component with*g* = *h* = 0 to be 1\*) indfd = Table[0, {*i*, 2\*Ng\*Nh}]; indfd[[Round[(Ng\*Nh + 1)/2]]] = 1.0; (\* forward incidence\*)

(\* Calculate the matrix Winv\*) exp*κ* = DiagonalMatrix[Table[Exp[*I*\**κ*a[[*i*, *i*]]\*dz\*2*π*/*λ*], {*i*, 4\*Ng\*Nh}]]; Winv = Dot[Tiinv, Ta, exp*κ*, Tainv, Ti]; (\* the total transfer matrix \*) *W* = Inverse[Winv];

(\* Calculate the diffraction and reflection fields \*) diff1 = Table[0, {*i*, 2\*Ng\*Nh}]; ref1 = Table[0, {*i*, 2\*Ng\*Nh}]; (\*initialize\*) W1 = *W*[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]]; W3 = *W*[[(2\*Ng\*Nh + 1);;4\*Ng\*Nh, 1;;2\*Ng\*Nh]]; diff1 = Dot[Inverse[W1], indfd]; (\* diffraction fields \*) ref1 = Dot[W3, Inverse[W1], indfd]; (\* reflection fields \*)

(\*Print diffraction and reflection fields as well as the corresponding *g*, *h* orders\*) Print["TE mode with multi-reflections"]; For[*i* = 1, *i* ≤ Ng\*Nh, Print[gindx[[*i*]], ", ", hindx[[*i*]], ", ", Abs[diff1[[*i*]]], ", ", Abs[ref1[[*i*]]]];*i*++; ];

(\*Solution(2) : solve diffractions and reflections without multi − reflections \*) (\*Calculate the matrixes related to virtual layer with *n* = (ne + no)/2\*) *ξ*avg = DiagonalMatrix[Table[Sqrt[((ne + no)/2.0)∧2 − nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]]; *ξ*avginv = DiagonalMatrix[Table[1.0/*ξ*avg[[*i*, *i*]], {*i*, Ng\*Nh}]]; T11 = nyd; T12 = nxd; T13 = nyd; T14 = nxd; T21 = Dot[nyd, *ξ*avg]; T22 = nI∧2\*Dot[nxd, *ξ*avginv]; T23 = −Dot[nyd, *ξ*avg]; T24 = −nI∧2\*Dot[nxd, *ξ*avginv]; T31 = −nxd; T32 = nyd; T33 = −nxd; T34 = nyd; T41 = Dot[nxd, *ξ*avg]; T42 = −nI∧2\*Dot[nyd, *ξ*avginv]; T43 = −Dot[nxd, *ξ*avg]; T44 = nI∧2\*Dot[nyd, *ξ*avginv]; T1i = Join[T11, T12, T13, T14, 2]; T2i = Join[T21, T22, T23, T24, 2]; T3i = Join[T31, T32, T33, T34, 2]; T4i = Join[T41, T42, T43, T44, 2]; Tavg = Join[T1i, T2i, T3i, T4i]; (\* transform matrix Ti\*) Tavginv = Inverse[Tavg]; (\* inverse of the transform matrix Ti \*) ClearAll[T11, T12, T13, T14, T21, T22, T23, T24, T31, T32, T33, T34, T41, T42, T43, T44]; ClearAll[T1i, T2i, T3i, T4i];

```
(* Calculate the transfer matrixes *)
S1 = Dot[Ta, expκ, Tainv];
Sent = Table[0, {i, 4*Ng*Nh}, {j, 4*Ng*Nh}];
Sext = Table[0, {i, 4*Ng*Nh}, {j, 4*Ng*Nh}];
W' = Inverse[Dot[Tavginv, Ti]];W1' = W'[[1;;2*Ng*Nh, 1;;2*Ng*Nh]];
W" = Inverse[Dot[Tiinv, Tavg]];W1" = W"[[1;;2*Ng*Nh, 1;;2*Ng*Nh]];
Sent[[1;;2*Ng*Nh, 1;;2*Ng*Nh]] = Inverse[W1']; Sent = Dot[Tavg, Sent];
Sext[[1;;2*Ng*Nh, 1;;2*Ng*Nh]] = Inverse[W1"]; Sext = Dot[Sext, Tavginv];
```
(\* Calculate the diffraction and reflection fields \*) indfd2 = Table[0, {*i*, 4\*Ng\*Nh}]; indfd2[[Round[(Ng\*Nh + 1)/2]]] = 1.0; (\* ignore backward \*) diff2 = Dot[Sext, S1, Sent, indfd2];

(\*Print diffraction and reflection fields as well as the corresponding *g*, *h* orders\*) Print["TE mode without multi-reflections"]; For[*i* = 1, *i* ≤ Ng\*Nh, Print[gindx[[*i*]], ", ", hindx[[*i*]], ", ", Abs[diff2[[*i*]]]];*i*++; ];

#### **5. References**

22 Will-be-set-by-IN-TECH

(\*set incident plane wave indfd,*e*.*g*. set the value of the component with*g* = *h* = 0 to be 1\*) indfd = Table[0, {*i*, 2\*Ng\*Nh}]; indfd[[Round[(Ng\*Nh + 1)/2]]] = 1.0; (\* forward incidence\*)

(\*Solution(1) : solve diffractions and reflections with multi − reflections\*)

exp*κ* = DiagonalMatrix[Table[Exp[*I*\**κ*a[[*i*, *i*]]\*dz\*2*π*/*λ*], {*i*, 4\*Ng\*Nh}]]; Winv = Dot[Tiinv, Ta, exp*κ*, Tainv, Ti]; (\* the total transfer matrix \*)

diff1 = Table[0, {*i*, 2\*Ng\*Nh}]; ref1 = Table[0, {*i*, 2\*Ng\*Nh}]; (\*initialize\*)

Print[gindx[[*i*]], ", ", hindx[[*i*]], ", ", Abs[diff1[[*i*]]], ", ", Abs[ref1[[*i*]]]];*i*++; ];

T31 = −nxd; T32 = nyd; T33 = −nxd; T34 = nyd; T41 = Dot[nxd, *ξ*avg];

(\*Print diffraction and reflection fields as well as the corresponding *g*, *h* orders\*)

(\*Solution(2) : solve diffractions and reflections without multi − reflections \*) (\*Calculate the matrixes related to virtual layer with *n* = (ne + no)/2\*)

T22 = nI∧2\*Dot[nxd, *ξ*avginv]; T23 = −Dot[nyd, *ξ*avg]; T24 = −nI∧2\*Dot[nxd, *ξ*avginv];

T42 = −nI∧2\*Dot[nyd, *ξ*avginv]; T43 = −Dot[nxd, *ξ*avg]; T44 = nI∧2\*Dot[nyd, *ξ*avginv];

ClearAll[T11, T12, T13, T14, T21, T22, T23, T24, T31, T32, T33, T34, T41, T42, T43, T44];

(\* Calculate the matrix Winv\*)

W1 = *W*[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]];

Print["TE mode with multi-reflections"];

− nx[[*i*, *i*]]∧2 − ny[[*i*, *i*]]∧2], {*i*, Ng\*Nh}]];

(\* Calculate the diffraction and reflection fields \*)

W3 = *W*[[(2\*Ng\*Nh + 1);;4\*Ng\*Nh, 1;;2\*Ng\*Nh]]; diff1 = Dot[Inverse[W1], indfd]; (\* diffraction fields \*) ref1 = Dot[W3, Inverse[W1], indfd]; (\* reflection fields \*)

*ξ*avg = DiagonalMatrix[Table[Sqrt[((ne + no)/2.0)∧2

Tavg = Join[T1i, T2i, T3i, T4i]; (\* transform matrix Ti\*)

*ξ*avginv = DiagonalMatrix[Table[1.0/*ξ*avg[[*i*, *i*]], {*i*, Ng\*Nh}]]; T11 = nyd; T12 = nxd; T13 = nyd; T14 = nxd; T21 = Dot[nyd, *ξ*avg];

T1i = Join[T11, T12, T13, T14, 2]; T2i = Join[T21, T22, T23, T24, 2]; T3i = Join[T31, T32, T33, T34, 2]; T4i = Join[T41, T42, T43, T44, 2];

Tavginv = Inverse[Tavg]; (\* inverse of the transform matrix Ti \*)

*W*' = Inverse[Dot[Tavginv, Ti]];W1' = *W*'[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]]; *W*" = Inverse[Dot[Tiinv, Tavg]];W1" = *W*"[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]]; Sent[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]] = Inverse[W1']; Sent = Dot[Tavg, Sent]; Sext[[1;;2\*Ng\*Nh, 1;;2\*Ng\*Nh]] = Inverse[W1"]; Sext = Dot[Sext, Tavginv];

*W* = Inverse[Winv];

For[*i* = 1, *i* ≤ Ng\*Nh,

ClearAll[T1i, T2i, T3i, T4i];

S1 = Dot[Ta, exp*κ*, Tainv];

(\* Calculate the transfer matrixes \*)

Sent = Table[0, {*i*, 4\*Ng\*Nh}, {*j*, 4\*Ng\*Nh}]; Sext = Table[0, {*i*, 4\*Ng\*Nh}, {*j*, 4\*Ng\*Nh}];


**5** 

*Algeria* 

**Wavelet Network Implementation on an** 

The approximation of general continuous functions by nonlinear networks is very useful for system modeling and identification. Such approximation methods can be used, for example, in black-box identification of nonlinear systems, signal processing, control, statistical data analysis, speech recognition, and artificial intelligence. Recently neural networks have been established as a general approximation tool for fitting nonlinear models from input/output data due to their ability of learning rather than complicated process functions (Gao, 2002). Their attractive property is the self-learning ability. A neural network can extract the system features from historical training data using the learning algorithm, requiring little or no a priori knowledge about the process (Patan, 2008). This is why during the past few years the nonlinear dynamic modelling of processes by neural networks has been extensively studied (Narendra & Parthasarathy, 1990; Nerrand et al., 1993; Levin, 1992; Rivals & Personnaz, 1996). In standard neural networks, the nonlinearities are approximated by superposition of

In the other hand, the wavelet theory has found many applications in function approximation, numerical analysis and signal processing. Though this attractive theory has offered efficient algorithms for various purposes, their implementations are usually limited to wavelets of small dimension. The reason is that constructing and storing wavelet basis of large dimension are of prohibitive cost. In order to handle problems of larger dimension, it is necessary to develop algorithms whose implementation is less sensitive to the dimension. And it is known that neural networks are powerful tools for handling problems of large

Due to the similarity between wavelet decomposition and one-hidden-layer neural networks, the idea of combining both wavelets and neural networks has been proposed in various works (Zhang & Benveniste, 1992; Pati & Krishnaprasad, 1993; Hong, 1992; Bakshi & Stephanopoulos, 1993; Tsatsanis & Giannakis, 1993;et al., 1994; Delyon et al., 1995; Saad Saoud & Khellaf, 2009). For example, in (Zhang & Benveniste, 1992) wavelet network is introduced as a class of feedforward networks composed of wavelets, in (Pati & Krishnaprasad, 1993) the discrete wavelet transform is used for analyzing and synthesizing

**1. Introduction** 

dimension.

sigmoidal functions (Cybenko, 1989).

**Inexpensive Eight Bit Microcontroller** 

Lyes Saad Saoud, Fayçal Rahmoune, Victor Tourtchine and Kamel Baddari

*Laboratory of Computer Science, Modeling, Optimization,* 

*Simulation and Electronic Systems (L.I.M.O.S.E), Department of Physics, Faculty of Sciences, University M'hamed Bougara Boumerdes* 


## **Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller**

Lyes Saad Saoud, Fayçal Rahmoune,

Victor Tourtchine and Kamel Baddari

*Laboratory of Computer Science, Modeling, Optimization, Simulation and Electronic Systems (L.I.M.O.S.E), Department of Physics, Faculty of Sciences, University M'hamed Bougara Boumerdes Algeria* 

### **1. Introduction**

24 Will-be-set-by-IN-TECH

86 Features of Liquid Crystal Display Materials and Processes

Lien, A.(1997). A detailed derivation of extended Jones matrix representation for twisted

Olivero, D & Oldano, C. (2003). Numerical methods for light propagation in large LC cells: a

Rokushima K. & Yamakita, J. (1983). Analysis of anisotropic dielectric gratings, *Journal of the*

Sutkowski M.; Grudniewski T.; Zmijan R.; Parka J. & Nowinowski K. E. (2006). Optical

Witzigmann, B; Regli, P & Fichtner, W. (1998). Rigorous electromagnetic simulation of liquid

Zhang, B. & Sheng, P. (2003). Optical measurement of azimuthal anchoring strength in nematic liquid crystals, *Physical Review E*, Vol. 67, Iss. 4, April 2003, pp. 041713-9.

new approach, *Liquid Crystals*, Vol. 30, Iss. 3, 2003, pp. 345-353.

*Optical Society of America*, Vol. 73, Iss. 7, July 1983, pp. 901-908.

171-175.

335-337.

pp.753-757.

nematic liquid crystal displays, *Liquid Crystals*, Vol. 22, No. 2, February 1997, pp.

data storage in LC cells, *Opto-Electronics Review*, Vol. 14, No. 4, December 2006, pp.

crystal displays, *Journal of the Optical Society of America A*, Vol. 15, Iss. 3, March 1998,

The approximation of general continuous functions by nonlinear networks is very useful for system modeling and identification. Such approximation methods can be used, for example, in black-box identification of nonlinear systems, signal processing, control, statistical data analysis, speech recognition, and artificial intelligence. Recently neural networks have been established as a general approximation tool for fitting nonlinear models from input/output data due to their ability of learning rather than complicated process functions (Gao, 2002). Their attractive property is the self-learning ability. A neural network can extract the system features from historical training data using the learning algorithm, requiring little or no a priori knowledge about the process (Patan, 2008). This is why during the past few years the nonlinear dynamic modelling of processes by neural networks has been extensively studied (Narendra & Parthasarathy, 1990; Nerrand et al., 1993; Levin, 1992; Rivals & Personnaz, 1996). In standard neural networks, the nonlinearities are approximated by superposition of sigmoidal functions (Cybenko, 1989).

In the other hand, the wavelet theory has found many applications in function approximation, numerical analysis and signal processing. Though this attractive theory has offered efficient algorithms for various purposes, their implementations are usually limited to wavelets of small dimension. The reason is that constructing and storing wavelet basis of large dimension are of prohibitive cost. In order to handle problems of larger dimension, it is necessary to develop algorithms whose implementation is less sensitive to the dimension. And it is known that neural networks are powerful tools for handling problems of large dimension.

Due to the similarity between wavelet decomposition and one-hidden-layer neural networks, the idea of combining both wavelets and neural networks has been proposed in various works (Zhang & Benveniste, 1992; Pati & Krishnaprasad, 1993; Hong, 1992; Bakshi & Stephanopoulos, 1993; Tsatsanis & Giannakis, 1993;et al., 1994; Delyon et al., 1995; Saad Saoud & Khellaf, 2009). For example, in (Zhang & Benveniste, 1992) wavelet network is introduced as a class of feedforward networks composed of wavelets, in (Pati & Krishnaprasad, 1993) the discrete wavelet transform is used for analyzing and synthesizing

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 89

The key problem in system identification is to find a suitable model structure within which a good model is to be found (Sjöberg et al. 1995). In general, the nonlinear dynamic system identification is the operation to determine a transformation operator *Ti*, for some desired

With *Tp(u)* and *Ti(u)* denote the system to be identified, which maps the compact set *U ϵ Rn*

The purpose is to find a class *Ti* such that *Tp* is represented by *Ti* adequately well. The operator *Tp* is defined by specific input-output pairs that are obtained from the inputs and the outputs of the system to be identified. Fig. 1 summarizes a system identification

Several nonlinear identification models can be found in the literature. The wavelet networks prove their capabilities (Zhang & Benveniste, 1992). The most important aspect of a wavelet network based identification scheme is the determination of an adaptive algorithm that minimizes the difference between the actual plant and the outputs of the identified model by using a set of training pairs which represent the approximate behavior of the actual

As mentioned in the introduction, it exists a great similarity between the wavelet network and the one hidden layer neural networks which gives it the opportunity to be an attractive alternative to the feed forward neural network or the Radial Basis Function (RBF) networks (Khalaf, 2006). The use of the wavelet as an activate function gives power to the network because it reduces naturally the noise and hence the resulting network becomes more simple due to the linear output. In this work, another aspect is given to this powerful architecture, which is the implementation of the whole wavelet network with its backpropagation

Wavelet function is a waveform that has a limited duration and an average value of zero. There exist several types of wavelet functions. In this work the Mexican Hat wavelet

1 4 2 2 ( ) 2 3 (1 )exp 2 (2)

(Abiyev, 2003) given by the equation 2 is used as a transfer function.

and *Y ϵ Rm*, and the identification model outputs respectively for the same input *u*.

() () , *Tu T u u U i p* (1)

**2. Nonlinear dynamic system identification** 

*ε>0*, so that (Efe & Kaynak, 1997):

Fig. 1. Wavelet networks identification

algorithm in the PIC microcontroller.

structure.

plant.

With:

feedforward neural networks, in (Hong, 1992) orthogonal wavelet bases are used for constructing wavelet-based neural network, and in (Saad Saoud & Khellaf, 2009) the dynamic wavelet networks is proposed and used to control the chemical reactor. Combining wavelets and neural networks can hopefully remedy the weakness of each other, resulting in networks with efficient constructive methods and capable of handling problems of moderately larger dimensions.

Hence, we can say that the neural network and the wavelet network are capable of modeling non-linear systems. On the basis of supplied training data the neural or the wavelet networks learn the relationship between the process input and output. The data have to be examined carefully before they can be used as a training set for network methods. The training sets consist of one or more input data and one or more output data (Roffel & Betlem, 2006). After the training of the network, a test-set of data should be used to verify whether the desired relationship was learned. These two operations (train and validate the network) are achieved generally by using the computer, the finding weights and bias are implemented either in the computer itself or through the implementation of the optimal network's parameters in the microcontroller (Gulbag et al., 2009; Cotton et al., 2008; Liung et al., 2003; Neelamegamand & Rajendran , 2005). One common drawback is that in both cases, in order to find the network's parameters precise calculations that are very processor intensive are required. This robust processing equipment can be expensive and rather large.

In several applications such as adaptive control (Plett, 2003), or predictive control (Liu et al., 1998), we need to adapt the network's parameters in real time and in this case the computer is very important to adapt the parameters. These problems were overcome with the implementation of the whole neural network with its backpropagation algorithm in the microcontroller, which is proposed in our previous work (Saad Saoud & Khellaf, 2011). In this later work a multilayer neural network is trained and validated using a very inexpensive and low end microcontroller, but the problems of larger dimensions still exist. For this case the real implementation of the whole wavelet network into an inexpensive microcontroller is proposed in this study.

The low end and inexpensive microcontroller PIC16F877A of Microchip trains and validates the wavelet network, and the well-known backpropagation algorithm is implemented to obtain the optimal network parameters. All the operations done by the microcontroller are shown through an alphanumeric liquid crystal display and several buttons are added in the embedded system which produces an ergonomic communication interface human/machine. The wavelet network takes more program memory place, for this reason the assembly language is preferred. The Continuous Stirred Tank Reactor (CSTR) system is chosen as a realistic nonlinear system to demonstrate the feasibility and the performance of the results found using the microcontroller. Several results will be presented in this chapter to give the reader more information about this field. A comparative study is made between the microcontroller and the computer.

The chapter is organized as follows: After the description of the nonlinear dynamic system identification in general and by using the wavelet network in particular, the implementation of the backpropagation algorithm for the wavelet network in the microcontroller. A comparison between wavelet network based on the eight bit microcontroller and those based on the computer is presented. To illustrate how effectively the eight bit microcontroller can learn nonlinear dynamic models, results for a Continuous Stirred Tank Reactor are given. All the electronic tools, electrical schemes and the implemented algorithm are discussed. The chapter concludes with few final remarks.

## **2. Nonlinear dynamic system identification**

88 Features of Liquid Crystal Display Materials and Processes

feedforward neural networks, in (Hong, 1992) orthogonal wavelet bases are used for constructing wavelet-based neural network, and in (Saad Saoud & Khellaf, 2009) the dynamic wavelet networks is proposed and used to control the chemical reactor. Combining wavelets and neural networks can hopefully remedy the weakness of each other, resulting in networks with efficient constructive methods and capable of handling problems of

Hence, we can say that the neural network and the wavelet network are capable of modeling non-linear systems. On the basis of supplied training data the neural or the wavelet networks learn the relationship between the process input and output. The data have to be examined carefully before they can be used as a training set for network methods. The training sets consist of one or more input data and one or more output data (Roffel & Betlem, 2006). After the training of the network, a test-set of data should be used to verify whether the desired relationship was learned. These two operations (train and validate the network) are achieved generally by using the computer, the finding weights and bias are implemented either in the computer itself or through the implementation of the optimal network's parameters in the microcontroller (Gulbag et al., 2009; Cotton et al., 2008; Liung et al., 2003; Neelamegamand & Rajendran , 2005). One common drawback is that in both cases, in order to find the network's parameters precise calculations that are very processor intensive are required. This robust processing equipment can be expensive and rather large. In several applications such as adaptive control (Plett, 2003), or predictive control (Liu et al., 1998), we need to adapt the network's parameters in real time and in this case the computer is very important to adapt the parameters. These problems were overcome with the implementation of the whole neural network with its backpropagation algorithm in the microcontroller, which is proposed in our previous work (Saad Saoud & Khellaf, 2011). In this later work a multilayer neural network is trained and validated using a very inexpensive and low end microcontroller, but the problems of larger dimensions still exist. For this case the real implementation of the whole wavelet network into an inexpensive

The low end and inexpensive microcontroller PIC16F877A of Microchip trains and validates the wavelet network, and the well-known backpropagation algorithm is implemented to obtain the optimal network parameters. All the operations done by the microcontroller are shown through an alphanumeric liquid crystal display and several buttons are added in the embedded system which produces an ergonomic communication interface human/machine. The wavelet network takes more program memory place, for this reason the assembly language is preferred. The Continuous Stirred Tank Reactor (CSTR) system is chosen as a realistic nonlinear system to demonstrate the feasibility and the performance of the results found using the microcontroller. Several results will be presented in this chapter to give the reader more information about this field. A comparative study is made between

The chapter is organized as follows: After the description of the nonlinear dynamic system identification in general and by using the wavelet network in particular, the implementation of the backpropagation algorithm for the wavelet network in the microcontroller. A comparison between wavelet network based on the eight bit microcontroller and those based on the computer is presented. To illustrate how effectively the eight bit microcontroller can learn nonlinear dynamic models, results for a Continuous Stirred Tank Reactor are given. All the electronic tools, electrical schemes and the implemented algorithm

moderately larger dimensions.

microcontroller is proposed in this study.

the microcontroller and the computer.

are discussed. The chapter concludes with few final remarks.

The key problem in system identification is to find a suitable model structure within which a good model is to be found (Sjöberg et al. 1995). In general, the nonlinear dynamic system identification is the operation to determine a transformation operator *Ti*, for some desired *ε>0*, so that (Efe & Kaynak, 1997):

$$\left\| T\_i(\mu) - T\_p(\mu) \right\| \le \varepsilon \quad \text{,} \mu \in \mathcal{U} \tag{1}$$

With *Tp(u)* and *Ti(u)* denote the system to be identified, which maps the compact set *U ϵ Rn* and *Y ϵ Rm*, and the identification model outputs respectively for the same input *u*.

The purpose is to find a class *Ti* such that *Tp* is represented by *Ti* adequately well. The operator *Tp* is defined by specific input-output pairs that are obtained from the inputs and the outputs of the system to be identified. Fig. 1 summarizes a system identification structure.

Several nonlinear identification models can be found in the literature. The wavelet networks prove their capabilities (Zhang & Benveniste, 1992). The most important aspect of a wavelet network based identification scheme is the determination of an adaptive algorithm that minimizes the difference between the actual plant and the outputs of the identified model by using a set of training pairs which represent the approximate behavior of the actual plant.

Fig. 1. Wavelet networks identification

As mentioned in the introduction, it exists a great similarity between the wavelet network and the one hidden layer neural networks which gives it the opportunity to be an attractive alternative to the feed forward neural network or the Radial Basis Function (RBF) networks (Khalaf, 2006). The use of the wavelet as an activate function gives power to the network because it reduces naturally the noise and hence the resulting network becomes more simple due to the linear output. In this work, another aspect is given to this powerful architecture, which is the implementation of the whole wavelet network with its backpropagation algorithm in the PIC microcontroller.

Wavelet function is a waveform that has a limited duration and an average value of zero. There exist several types of wavelet functions. In this work the Mexican Hat wavelet (Abiyev, 2003) given by the equation 2 is used as a transfer function.

$$\Psi(\pi) = \begin{pmatrix} 2 \end{pmatrix} \begin{pmatrix} \sqrt{3} \\ \end{pmatrix} \text{ } \pi^{-1/4} (1 - \pi^2) \exp\left(-\pi^2 / 2\right) \tag{2}$$

With:

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 91

**Dilation d=0.5**

t=0

t=1

(4)

*ek kk y y*ˆ (5)

*ww e new old <sup>k</sup>* (6)

*new old <sup>k</sup> t t w he* (7)

*new old <sup>k</sup> d d w he* (8)

1 4 2 2 2 3 (3 )exp 2 *old h d* (9)


**Input (u(k))**

2 1

*T*

*T*

Fig. 3. Mexican Hat wavelet function for different translations and a dilation equals to 0.5

1 2 *m k k E e* 

Assuming that the objective is to minimize this kind of performance function, the network

As we can see from the above equations, the updating network parameters are performed by evaluating the gradient of the performance function with respect to each individual

parameters updating rule is given by the equations (6) to (9):


network parameter in the network.

With

With:


0

0.2

**Mexican Hat wavelet function** 

**()**

0.4

0.6

0.8

t=-1

1

$$\tau = \frac{u\_i(k) - t\_i}{d\_i}$$

The network output is simply the sum of the neuron's outputs of the hidden layer and it is given by the equation (3) as follows:

$$
\hat{\boldsymbol{y}}\left(\boldsymbol{k}\right) = \sum\_{i=1}^{n} w\_{i} \boldsymbol{\upmu}\left(\boldsymbol{\uptau}\right) \tag{3}
$$

Where :

*u(k)* is the neuron input vector of the dimension [*nx1*].

*y(k)* is the neuron output at time instant *k*.

*ψ* is wavelet activation function of the neuron with a translation *t* and a dilation *d*.

Figures 2 and 3 are the representation of the Mexican Hat wavelet function for different translations and dilations:

Fig. 2. Mexican Hat wavelet function for different dilations and a translation equals to zero

The backpropagation algorithm is the most popular training method which is widely used in the neural network applications (Efe, 1996). The backpropagation algorithm is, however, very general and not limited to one-hidden-layer sigmoid neural network models. Instead, it could be applied to all network models (Sjöberg et al., 1995).

In this chapter we use this elegant technique to train the wavelet network. The method minimizes the performance or the so-called the cost function defined on the actual and desired outputs of the network by the equation (4). When updating each individual network parameter (weight, translation and dilation), the gradient information obtained from the differentiation of the cost function is used. As a matter of fact, we are looking for least mean squares. This can be attained by moving the network parameters vectors in a direction such that the performance function decreases. It is obvious that this direction is the negative gradient direction of the performance function.

Fig. 3. Mexican Hat wavelet function for different translations and a dilation equals to 0.5

$$E = \frac{1}{2} \sum\_{k=1}^{m} e\_k^2 \tag{4}$$

With

90 Features of Liquid Crystal Display Materials and Processes

The network output is simply the sum of the neuron's outputs of the hidden layer and it is

 1

Figures 2 and 3 are the representation of the Mexican Hat wavelet function for different

**Translation t=0**

d=0.5 d=1

d=2

d=0.1


**Input (u(k))**

Fig. 2. Mexican Hat wavelet function for different dilations and a translation equals to zero The backpropagation algorithm is the most popular training method which is widely used in the neural network applications (Efe, 1996). The backpropagation algorithm is, however, very general and not limited to one-hidden-layer sigmoid neural network models. Instead, it

In this chapter we use this elegant technique to train the wavelet network. The method minimizes the performance or the so-called the cost function defined on the actual and desired outputs of the network by the equation (4). When updating each individual network parameter (weight, translation and dilation), the gradient information obtained from the differentiation of the cost function is used. As a matter of fact, we are looking for least mean squares. This can be attained by moving the network parameters vectors in a direction such that the performance function decreases. It is obvious that this direction is the negative

could be applied to all network models (Sjöberg et al., 1995).

gradient direction of the performance function.

*i i yk w* 

(3)

<sup>ˆ</sup> *<sup>n</sup>*

*ψ* is wavelet activation function of the neuron with a translation *t* and a dilation *d*.

given by the equation (3) as follows:

*y(k)* is the neuron output at time instant *k*.



0

0.2

**Mexican Hat wavelet function** 

**()**

0.4

0.6

0.8

1

translations and dilations:

*u(k)* is the neuron input vector of the dimension [*nx1*].

Where :

( ) *i i i uk t d* 

$$e\_k = \left(\mathcal{y}\_k - \hat{\mathcal{y}}\_k\right) \tag{5}$$

Assuming that the objective is to minimize this kind of performance function, the network parameters updating rule is given by the equations (6) to (9):

$$
\omega\_{new} = \omega\_{old} + \eta \,\mu\left(\tau\right) e\_k \tag{6}
$$

$$t\_{new} = t\_{old} + \eta \, w^T h \, e\_k \tag{7}$$

$$d\_{new} = d\_{old} + \eta \,\mathrm{tr} \, w^T h \, e\_k \tag{8}$$

With:

$$\hbar h = \left(2\sqrt{\sqrt{3}}\right) \pi^{-1/4} \pi (3 - \pi^2) \exp\left(-\pi^2 \left/2\right)\right) d\_{old} \tag{9}$$

As we can see from the above equations, the updating network parameters are performed by evaluating the gradient of the performance function with respect to each individual network parameter in the network.

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 93

;\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\* ; ; This is a sample program of the wavelet network approach written in assembly language. ; ; © Lyes Saad Saoud, Fayçal Rahmoune, Victor Tourtchine & Kamel Baddari, 2011 ; ;\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\* ;

;-------------------------------------------------------------------------------------------------------------------------; ; For read 500 values from the extern EEPROM memory ; ;-------------------------------------------------------------------------------------------------------------------------;

 CALL LOAD\_DATA ; Load one input, output regressor CALL WAVENET\_OUTPUT ; Compute the wavelet network output CALL ERROR ; Compute the error between the ; modeled and estimated outputs CALL ADAPTATION ; Call the backpropagation subprogram CALL MSE ; Compute the last mean squared error

 DECFSZ Counter1,f ; Decrease the counter 1 ; If counter 1 equal zero jump GOTO Loop1 ; Else, loop to counter 1 DECFSZ Counter2,f ; Decrease counter 2

BCF PORTC,1 ; Orient LCD for control

BSF PORTC, 1 ; Orient LCD for data

 SENDCAR 'S' ; for display SENDCAR 'E' ; MSE= SENDCAR '=' ;

SENDCAR 'M' ; Send the letters "M, S, E, ="

 ; If counter 2 equal zero jumps GOTO Loop2 ; Else, loop to counter 2

;----------------------------------------------------------------------------------------------------------------------;

; instruction at the end

 SENDCAR B'11000000' ; SENDCAR is a macro to send data ; to the display, this data tell the ; display to start writing in the line 2 ; This macro always validate the

; procedure

 MOVLW .100 ; 100 iterations for the ; backpropagation

MOVWF Counter3 ; In Counter3

 MOVLW .5 ; For 5 times MOVWF Counter2 ; In Counter2

 MOVLW .100 ; For 100 times MOVWF Counter1 ; In Counter1

START . . .

Loop3

Loop2

Loop1

## **3. PIC Implementation of the full wavelet network**

In our previous paper (Saad Saoud & Khellaf, 2011) the neural network is implemented and trained using Microchip microcontroller (PIC16F876A). In this work, the extension of this card is used to train and validate another more powerful model which is the wavelet network. In this work, the microcontroller type is replaced by the PIC16F877A to give more power and utility to the final realized card.

First of all, why did we choose this type of microcontroller?

The PIC microcontrollers not only are the most widely used and well known microcontrollers, they are also the best supported. In fact, PIC system design and programming has become a powerful specialization with a large number of professional and amateur specialists. There are hundreds of websites devoted to PIC-related topics. An entire cottage industry of PIC software and hardware has flourished around this technology (Sanchez & Canton, 2007).

In this chapter and like our previous work we keep the mid-range PIC microcontroller, but we move to an upgraded version: the PIC16F877A, to give the new designed card better extensions, such as, more inputs outputs giving the user more flexibility when using the card in the control or other real practical situations.

The wavelet network implementation is written in assembly language based on the MPLAB tool distributed freely by Microchip Company. A sample of the 9162 assembly lines taken as the principal program is shown in the figure 4. The wavelet network hexadecimal source (Wavelet.hex) can be found in Intech publisher. Like the old realized card, the microcontroller has to find the optimal wavelet network parameters by using the well known backpropagation algorithm based on the memorized input output pairs in the EEPROM memory. The photo of the realized card based on the PIC16F877A microcontroller is shown in figure 5.

As shown in the program sample, the loops (Loop 1 and Loop 2) are used to read the 500 pairs input-output. The "Loop 3" in Fig. 4 is simply used to train the network. It can be noticed that the program is very flexible, and it gives the designers the opportunity to change the iteration number, the number of the pairs of modeling data and also the type of the data itself, just they have to respect the input regressor.

#### **4. Practical results**

To test the performance of the realized card (Fig.5), the Continuous Stirred Tank Reactor (CSTR) chemical reactor is chosen. The secrets of this system are the high nonlinearity and simplicity of the mathematical model. For this reason several researchers use the CSTR as a simulated or practical plant to validate their results such as the references (Lightbody & Irwin , 1997; Morningred et al.,1990; Henson & Seborg , 1990; Espinosa et al., 2005; Saad Saoud & Khellaf, 2009; Saad Saoud et al., 2011), and it consists simply of an irreversible, exothermic reaction. A→B, in a constant volume rate cooled by a single coolant stream which can be modeled by the following equations :

$$
\dot{\mathbf{C}}\_a(t) = \frac{q}{v} (\mathbf{C}\_{a\nu} - \mathbf{C}\_a(t)) - k\_0 \mathbf{C}\_a(t) e^{-\frac{E}{RT(t)}} \tag{10}
$$

$$\dot{T}(t) = \frac{q}{v}(T\_0 - T(t)) - k\_1 \mathbf{C}\_a(t) e^{-\frac{E}{RT(t)}} + k\_2 q\_c(t) (1 - e^{-\frac{k\_3}{q\_c(t)}}) (T\_{c0} - T(t)) \tag{11}$$

In our previous paper (Saad Saoud & Khellaf, 2011) the neural network is implemented and trained using Microchip microcontroller (PIC16F876A). In this work, the extension of this card is used to train and validate another more powerful model which is the wavelet network. In this work, the microcontroller type is replaced by the PIC16F877A to give more

The PIC microcontrollers not only are the most widely used and well known microcontrollers, they are also the best supported. In fact, PIC system design and programming has become a powerful specialization with a large number of professional and amateur specialists. There are hundreds of websites devoted to PIC-related topics. An entire cottage industry of PIC software and hardware has flourished around this technology

In this chapter and like our previous work we keep the mid-range PIC microcontroller, but we move to an upgraded version: the PIC16F877A, to give the new designed card better extensions, such as, more inputs outputs giving the user more flexibility when using the

The wavelet network implementation is written in assembly language based on the MPLAB tool distributed freely by Microchip Company. A sample of the 9162 assembly lines taken as the principal program is shown in the figure 4. The wavelet network hexadecimal source (Wavelet.hex) can be found in Intech publisher. Like the old realized card, the microcontroller has to find the optimal wavelet network parameters by using the well known backpropagation algorithm based on the memorized input output pairs in the EEPROM memory. The photo of

As shown in the program sample, the loops (Loop 1 and Loop 2) are used to read the 500 pairs input-output. The "Loop 3" in Fig. 4 is simply used to train the network. It can be noticed that the program is very flexible, and it gives the designers the opportunity to change the iteration number, the number of the pairs of modeling data and also the type of

To test the performance of the realized card (Fig.5), the Continuous Stirred Tank Reactor (CSTR) chemical reactor is chosen. The secrets of this system are the high nonlinearity and simplicity of the mathematical model. For this reason several researchers use the CSTR as a simulated or practical plant to validate their results such as the references (Lightbody & Irwin , 1997; Morningred et al.,1990; Henson & Seborg , 1990; Espinosa et al., 2005; Saad Saoud & Khellaf, 2009; Saad Saoud et al., 2011), and it consists simply of an irreversible, exothermic reaction. A→B, in a constant volume rate cooled by a single coolant stream

<sup>0</sup> ( ) ( ( )) ( )

*<sup>q</sup> RT t C t C C t kC te <sup>a</sup> <sup>v</sup> ao a <sup>a</sup>*

01 2 0 ( ) ( ( )) ( ) ( )(1 )( ( )) *<sup>c</sup>*

*<sup>q</sup> RT t q t Tt T Tt kC te k <sup>v</sup> ac c <sup>q</sup> t e T Tt*

( )

3

*E*

 (10)

( ) ( )

*E k*

(11)

the realized card based on the PIC16F877A microcontroller is shown in figure 5.

**3. PIC Implementation of the full wavelet network** 

First of all, why did we choose this type of microcontroller?

card in the control or other real practical situations.

the data itself, just they have to respect the input regressor.

which can be modeled by the following equations :

power and utility to the final realized card.

(Sanchez & Canton, 2007).

**4. Practical results** 

;\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\* ; ; This is a sample program of the wavelet network approach written in assembly language. ; ; © Lyes Saad Saoud, Fayçal Rahmoune, Victor Tourtchine & Kamel Baddari, 2011 ; ;\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\* ; START . . . MOVLW .100 ; 100 iterations for the ; backpropagation ; procedure MOVWF Counter3 ; In Counter3 Loop3 ;-------------------------------------------------------------------------------------------------------------------------; ; For read 500 values from the extern EEPROM memory ; ;-------------------------------------------------------------------------------------------------------------------------; MOVLW .5 ; For 5 times MOVWF Counter2 ; In Counter2 Loop2 MOVLW .100 ; For 100 times MOVWF Counter1 ; In Counter1 Loop1 CALL LOAD\_DATA ; Load one input, output regressor CALL WAVENET\_OUTPUT ; Compute the wavelet network output CALL ERROR ; Compute the error between the ; modeled and estimated outputs CALL ADAPTATION ; Call the backpropagation subprogram CALL MSE ; Compute the last mean squared error DECFSZ Counter1,f ; Decrease the counter 1 ; If counter 1 equal zero jump GOTO Loop1 ; Else, loop to counter 1 DECFSZ Counter2,f ; Decrease counter 2 ; If counter 2 equal zero jumps GOTO Loop2 ; Else, loop to counter 2 ;----------------------------------------------------------------------------------------------------------------------; BCF PORTC,1 ; Orient LCD for control SENDCAR B'11000000' ; SENDCAR is a macro to send data ; to the display, this data tell the ; display to start writing in the line 2 ; This macro always validate the ; instruction at the end BSF PORTC, 1 ; Orient LCD for data SENDCAR 'M' ; Send the letters "M, S, E, =" SENDCAR 'S' ; for display SENDCAR 'E' ; MSE= SENDCAR '=' ;

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 95

Fig. 5. The realized card based on the PIC microcontroller using the wavelet network

These operations and the others are summarized in the following algorithm:

3. Calculate the wavelet network output using the equation (3)

Let's take for example the concentration *Ca=0.11944 mol/l*.

4. Calculate the error using the equation (5)

1. Initialize the weights and translations randomly small and the dilations by ones.

5. Update the weights, the translations and the dilations with the equations (6) to (9). 6. If the *E < ε* or the number of the iteration is achieved, End the operation

We have to normalize this value between 0.1 and 0.9, for this we make this operation:

To illustrate this algorithm and enter deeply into the microcontroller program memory we

denominators.

2. Load the data

Else goto the step 2.

give the following example :

The data used for modeling can be found in (De Moor, 1998) and it is shown in figure 7. We choose randomly a sequence of 750 samples of the input-output from the whole given data. As we have made in our previous works (Saad Saoud & Khellaf, 2011; Saad Saoud et al., 2011), before using the data, we have to make several changes to be accepted by the network. The used data saved in the EEPROM passed through two necessary operations, first it has been normalized to 0.1 and 0.9 and second it is multiplied by 36408 to cover the maximum range and not exceeding the positive signed range of 16 bit (7FFF hex) of data. It should be mentioned that in very few cases, error values could be displayed on the LCD and this is due the floating point format difference between the computer and the PIC microcontroller (PIC 16F877A, 24 bit floating point format). In the microcontroller, all these operations will be reversed to have real values. Also it can be seen in the previous sample assembly language program that 500 samples were used to train the network and the rest of data is used for the network validation. At the beginning, the microcontroller trains the network to find the optimal parameters. All the parameters except the dilations were initialized randomly, hence the network dilations were initialized ones to avoid zeros in the



The process describes the reaction of two products, which are mixed and react to generate a compound *A* having a concentration *Ca(t)*, with the temperature of the mixture *T(t)*. This reaction is exothermic. The generated heat acts to slow the reaction. The reaction is controlled by introducing a coolant flow rate *qc(t)*, which helps to change the temperature and thereby the concentration. *Ca0* is the inlet feed concentration, *q* is the process flow rate, *T0* and *Tc0* are the inlet feed and coolant temperatures. All these values are assumed constant at nominal values. In the same way, *k0, E/R, υ, k1, k2* and *k3* are thermodynamic and chemical constants.

$$k\_1 = -\frac{\Delta H k\_0}{\rho \mathcal{C}\_p} \qquad \qquad k\_2 = -\frac{\rho\_c \mathcal{C}\_{pc}}{\rho \mathcal{C}\_p \upsilon} \qquad \qquad k\_3 = -\frac{h\_a}{\rho\_c \mathcal{C}\_{pc}}.$$

The nominal conditions for a product concentration *Ca = 0.1 mol/l* are: *T = 438.54K, qc = 103.41 l/min.* 

.

Fig. 5. The realized card based on the PIC microcontroller using the wavelet network

The data used for modeling can be found in (De Moor, 1998) and it is shown in figure 7. We choose randomly a sequence of 750 samples of the input-output from the whole given data. As we have made in our previous works (Saad Saoud & Khellaf, 2011; Saad Saoud et al., 2011), before using the data, we have to make several changes to be accepted by the network. The used data saved in the EEPROM passed through two necessary operations, first it has been normalized to 0.1 and 0.9 and second it is multiplied by 36408 to cover the maximum range and not exceeding the positive signed range of 16 bit (7FFF hex) of data. It should be mentioned that in very few cases, error values could be displayed on the LCD and this is due the floating point format difference between the computer and the PIC microcontroller (PIC 16F877A, 24 bit floating point format). In the microcontroller, all these operations will be reversed to have real values. Also it can be seen in the previous sample assembly language program that 500 samples were used to train the network and the rest of data is used for the network validation. At the beginning, the microcontroller trains the network to find the optimal parameters. All the parameters except the dilations were initialized randomly, hence the network dilations were initialized ones to avoid zeros in the denominators.

These operations and the others are summarized in the following algorithm:


94 Features of Liquid Crystal Display Materials and Processes

CALL Display\_MSE ; Calculate and display the mean

MOVLW 0x0E ; Point to the first output address

MOVLW 0x14 ; Point to the next output address

;----------------------------------------------------------------------------------------------------------------------; ; Validate with new input-output pairs and read the trained and validated data from ; ; the EEPROM memory and display them through the alphanumerical LCD ; ;----------------------------------------------------------------------------------------------------------------------;

; The initial input address

 DECFSZ Counter3,f ; Decrease the counter 3 ; If counter 3 equal zero jumps

GOTO Loop3 ; Else, loop to counter 3 . .

 CALL Validation\_Part ; Validate the founding model CALL Read\_data ; Read the trained and validated ; data from the EEPROM memory ; and display them through the

; alphanumerical LCD

Fig. 4. A sample program of the assembly wavelet network implementation

GOTO START ; Start from the beginning . .

In the same way, *k0, E/R, υ, k1, k2* and *k3* are thermodynamic and chemical constants.

*k*

0

The nominal conditions for a product concentration *Ca = 0.1 mol/l* are:

*p*

1

*T = 438.54K, qc = 103.41 l/min.* 

*Hk <sup>k</sup> C* <sup>2</sup>

The process describes the reaction of two products, which are mixed and react to generate a compound *A* having a concentration *Ca(t)*, with the temperature of the mixture *T(t)*. This reaction is exothermic. The generated heat acts to slow the reaction. The reaction is controlled by introducing a coolant flow rate *qc(t)*, which helps to change the temperature and thereby the concentration. *Ca0* is the inlet feed concentration, *q* is the process flow rate, *T0* and *Tc0* are the inlet feed and coolant temperatures. All these values are assumed constant at nominal values.

> *c pc p C*

<sup>3</sup> *<sup>a</sup>*

*c pc*

*<sup>h</sup> <sup>k</sup> <sup>C</sup>*

*C*

; equal to 0x0E14

 ; squared error MSE CLRF MSEEXP ; Clear register MSEEXP CLRF MSEAARGB0 ; Clear register MSEAARGB0 CLRF MSEAARGB1 ; Clear register MSEAARGB1 CLRF eepm\_output ; Point to the first output address CLRF eepm\_output+1 ; Point to the next output address ; The initial output address equal to

; zero

MOVWF eepm\_input ;

MOVWF eepm\_input+1 ;

.

.


To illustrate this algorithm and enter deeply into the microcontroller program memory we give the following example :

Let's take for example the concentration *Ca=0.11944 mol/l*.

We have to normalize this value between 0.1 and 0.9, for this we make this operation:

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 97

0 100 200 300 400 500 600 700

0 100 200 300 400 500 600 700

Fig. 7. Data used to train and validate the wavelet network using the PIC microcontroller

N° of samples

Fig. 6. The electronic circuit of the realized card

90

0.06 0.08 0.1 0.12 0.14 0.16

Concentration (mol/l)

95

100

Flow rate (l/h)

105

110

*Ca\_Normalized= a\_NxCa + b\_N* 

When the coefficients *a\_N* and *b\_N* are the normalization coefficients for the whole concentration data and in this case are :

*a\_N=9.2937 and b\_N=-0.4638* 

We can find simply *Ca\_Normalized= 0.64625*.

And before saving this value in the EEPROM we have to multiply it by 36408. So the value that should be saved is 5BE8 hex. This operation will be done with all the taken samples.

The microcontroller will make the inverse operations at each time it reads the values from the EEPROM.

In the realized card shown in the figure 5, we can find, in addition to the microcontroller, the EEPROM 24C512, the power supply part and the alphanumeric Liquid Cristal Display (LCD). This last one has a capital importance, it is the bridge between the users and the embedded card and without it the user cannot manipulate the card. Its use gives the realized card power and simplicity. It displays all the necessary operations. For example, the user can start and see the training part of the network, the mean squared error will be displayed at each time the program finishes the three loops (Loop 1 to 3) in the sample program of the figure 4. The users can read the validated data memorized in the EEPROM at the end of the training part and at each time the system can be interrupted by pushing on the rest button.

As shown in the electronic circuit (Fig. 6) it contains few components, and we can say that the deserved work can be done with three principal components: the microcontroller with its basic circuits, the memory and the LCD.

The comparison between the obtained results using two instruments (the computer and the microcontroller) and based on the extensive practical results carried out within the course of this study is shown by figures 8 and 9 and it is given in the table 1. As remarked, the obtained results using the realized card is very close to the real results, and sometimes better than the results obtained using the computer especially in the validation part (Fig. 9). This is due to the well known problem of the wavelet network, which is the choice of the initial wavelet network parameters. These results give the microcontroller applications an expansion in the artificial intelligence field.

On the other hand, when we want to compare the card based on the neural network which has been realized in our previous work and the card presented in this chapter, we can find several points. First, the assembly language line number is smaller in the wavelet network based card than the program based on the neural network. This difference gives the card an important consequence which is the time reduction 2.77 second which is greater in the card based on the neural network (Saad Saoud & Khellaf, 2011). These two important advantages are the result of the network architecture itself, because the wavelet network has a linear output which makes the calculations easier.

Second, the new card gives the designer more simplicity of use, because the assembly program is written in the way that we can make a call of the standard subroutines, and we can also change at any time the network parameters such as the number of samples for training and validation or the number of iterations. We can say also that there are several disadvantages of the new card but not so important. The card cost is a little more expensive that the old card (Saad Saoud & Khellaf, 2011), and the wavelet network needs a special care during parameters initialization.

When the coefficients *a\_N* and *b\_N* are the normalization coefficients for the whole

And before saving this value in the EEPROM we have to multiply it by 36408. So the value that should be saved is 5BE8 hex. This operation will be done with all the taken samples. The microcontroller will make the inverse operations at each time it reads the values from

In the realized card shown in the figure 5, we can find, in addition to the microcontroller, the EEPROM 24C512, the power supply part and the alphanumeric Liquid Cristal Display (LCD). This last one has a capital importance, it is the bridge between the users and the embedded card and without it the user cannot manipulate the card. Its use gives the realized card power and simplicity. It displays all the necessary operations. For example, the user can start and see the training part of the network, the mean squared error will be displayed at each time the program finishes the three loops (Loop 1 to 3) in the sample program of the figure 4. The users can read the validated data memorized in the EEPROM at the end of the training part and at each time the system can be interrupted by pushing on

As shown in the electronic circuit (Fig. 6) it contains few components, and we can say that the deserved work can be done with three principal components: the microcontroller with

The comparison between the obtained results using two instruments (the computer and the microcontroller) and based on the extensive practical results carried out within the course of this study is shown by figures 8 and 9 and it is given in the table 1. As remarked, the obtained results using the realized card is very close to the real results, and sometimes better than the results obtained using the computer especially in the validation part (Fig. 9). This is due to the well known problem of the wavelet network, which is the choice of the initial wavelet network parameters. These results give the microcontroller applications an

On the other hand, when we want to compare the card based on the neural network which has been realized in our previous work and the card presented in this chapter, we can find several points. First, the assembly language line number is smaller in the wavelet network based card than the program based on the neural network. This difference gives the card an important consequence which is the time reduction 2.77 second which is greater in the card based on the neural network (Saad Saoud & Khellaf, 2011). These two important advantages are the result of the network architecture itself, because the wavelet network has a linear

Second, the new card gives the designer more simplicity of use, because the assembly program is written in the way that we can make a call of the standard subroutines, and we can also change at any time the network parameters such as the number of samples for training and validation or the number of iterations. We can say also that there are several disadvantages of the new card but not so important. The card cost is a little more expensive that the old card (Saad Saoud & Khellaf, 2011), and the wavelet network needs a special care

*Ca\_Normalized= a\_NxCa + b\_N* 

*a\_N=9.2937 and b\_N=-0.4638* 

the EEPROM.

the rest button.

concentration data and in this case are :

We can find simply *Ca\_Normalized= 0.64625*.

its basic circuits, the memory and the LCD.

expansion in the artificial intelligence field.

output which makes the calculations easier.

during parameters initialization.

Fig. 6. The electronic circuit of the realized card

Fig. 7. Data used to train and validate the wavelet network using the PIC microcontroller

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 99

Reference model Computer model Microcontroller model

500 600 700

N° of samples

Fig. 9. Validation results of the simulated CSTR reactor using neural network based on the



Hence, we can say, the choice of the LCD type is based on the practical situation in where

In this chapter, a direct practical application of the Liquid Cristal Display (LCD) is presented and proved by the real implementation of the full wavelet network with its backpropagation in a very inexpensive microcontroller. The realized cad is embedded and very simple to use by people having a small knowledge in the electronic and the assembly language programming. The embedded card is tested using the celebrated nonlinear system which is the CSTR chemical reactor. And in all cases, the programmer can at any time change the system to be modeled using the wavelet network respecting the input regressor and the data manipulation before saving it. In our future work, we try to use the realized card in real

In this appendix we try to give the reader of this chapter the tools to realize and test by

First, it is very important to realize the card, for this reason the card's PCB is in pdf format that can be found within this book (wavelet mask.pdf and wavelet card.pdf). Second, the

users should implement the components like it is shown in figures 5 and 6.

hand, bigger LCD occupies more place which means large card dimension.

control situations and giving the architecture a dynamic aspect.

0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14

display type.

**6. Conclusion** 

**7. Appendix** 

microcontroller and the computer.

the realized card will be implemented.

himself the proposed card in this chapter.

Cencentration (mol/l)

Fig. 8. Comparison results between computer and microcontroller for the training part of the CSTR using wavelet network based.

## **5. Moving to other LCD technologies**

In this chapter, a direct Liquid Crystal Display (LCD) application is illustrated. Even though it is used as a slave of the PIC16F877A, it could be noted that without this tool the electronic card can not be used. This proves the usefulness of the liquid crystal display in the industrial life.

This card is normally operated with any LCD type, which is actually in our work the Twisted Nematic (TN)-Display. But it can be improved with other display types such us the In-Plane Switching (IPS)- and Multidomain–Vertical-Alignment (MVA)-Displays. In this work, we take the reference (Willem den Boer, 2005), as an example to describe these modes.


Table 1. Comparison between the two strategies (computer and microcontroller)

An often-quoted drawback of conventional TN LCDs has been their poor viewing angle behavior. Several dramatic improvements in viewing angle have been developed over the past few years. The most important ones are (Willem den Boer, 2005): The In-Plane-Switching (IPS) LC mode and the Multidomain-Vertical-Alignment (MVA) LC mode.

In both cases, the aim is to obtain a good viewing angle. The use of one of them gives the realized card another aspect, and improves the quality of the electronic system. Like all systems, moving to another type has advantages and disadvantages. We can cite here the principal disadvantages:

Fig. 9. Validation results of the simulated CSTR reactor using neural network based on the microcontroller and the computer.


Hence, we can say, the choice of the LCD type is based on the practical situation in where the realized card will be implemented.

## **6. Conclusion**

98 Features of Liquid Crystal Display Materials and Processes

0 100 200 300 400 500

N° of samples

Fig. 8. Comparison results between computer and microcontroller for the training part of the

In this chapter, a direct Liquid Crystal Display (LCD) application is illustrated. Even though it is used as a slave of the PIC16F877A, it could be noted that without this tool the electronic card can not be used. This proves the usefulness of the liquid crystal display in the

This card is normally operated with any LCD type, which is actually in our work the Twisted Nematic (TN)-Display. But it can be improved with other display types such us the In-Plane Switching (IPS)- and Multidomain–Vertical-Alignment (MVA)-Displays. In this work, we take the reference (Willem den Boer, 2005), as an example to describe these modes.

> **Computer Pentium Dual-core inside speed 1.73 GHz**

Table 1. Comparison between the two strategies (computer and microcontroller)

An often-quoted drawback of conventional TN LCDs has been their poor viewing angle behavior. Several dramatic improvements in viewing angle have been developed over the past few years. The most important ones are (Willem den Boer, 2005): The In-Plane-Switching (IPS) LC mode and the Multidomain-Vertical-Alignment (MVA) LC mode. In both cases, the aim is to obtain a good viewing angle. The use of one of them gives the realized card another aspect, and improves the quality of the electronic system. Like all systems, moving to another type has advantages and disadvantages. We can cite here the

**Instrument's cost** Expensive Inexpensive **Space occupation** Big space consumed Embedded

**Wavelet network trained through** 

0.183 for 500 samples 2.77 for 500 samples

**Microcontroller PIC16F877A speed 20 MHz** 

Reference model Computer model Microcontroller model

0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12

CSTR using wavelet network based.

industrial life.

**One iteration execution time (sec)** 

principal disadvantages:

**5. Moving to other LCD technologies** 

Cencentration (mol/l)

In this chapter, a direct practical application of the Liquid Cristal Display (LCD) is presented and proved by the real implementation of the full wavelet network with its backpropagation in a very inexpensive microcontroller. The realized cad is embedded and very simple to use by people having a small knowledge in the electronic and the assembly language programming. The embedded card is tested using the celebrated nonlinear system which is the CSTR chemical reactor. And in all cases, the programmer can at any time change the system to be modeled using the wavelet network respecting the input regressor and the data manipulation before saving it. In our future work, we try to use the realized card in real control situations and giving the architecture a dynamic aspect.

## **7. Appendix**

In this appendix we try to give the reader of this chapter the tools to realize and test by himself the proposed card in this chapter.

First, it is very important to realize the card, for this reason the card's PCB is in pdf format that can be found within this book (wavelet mask.pdf and wavelet card.pdf). Second, the users should implement the components like it is shown in figures 5 and 6.

Wavelet Network Implementation on an Inexpensive Eight Bit Microcontroller 101

Kreinovich, V.; Sirisaengtaksin, O. & Cabrera, S. (1994). Wavelet neural networks are

Levin, A. U. (1992). Neural networks in dynamical systems: a system theoretic approach,

Lightbody, G. & Irwin, G. (1997). Nonlinear control structures based on embedded neural system models. *IEEE Transaction on Neural Networks*, 8, pp. 553–567. Liu, G. P.; Kadirkamanathan, V. & Billings, S. A. (1998). Predictive control for non-linear

Liung, T. K.; Mashor, M. Y.; Isa, N. A. M.; Ali, A. N. & Othman, N. H. (2003). Design of a

Morningred, J. D.; Paden, B. E.; Seborg D. E. & Mellichamp, D. A.. An Adaptive Nonlinear

Narendra, K. S. & Parthasarathy, K. (1990). Identification and Control of Dynamical Systems Using Neural Networks, *IEEE Transactions on Neural Networks*, Vol.1 (1), pp. 4-27. Neelamegamand, P. & Rajendran, A. (2005). Neural Network Based Density Measurement,

Nerrand, O.; Roussel-Ragot, P.; Personnaz, L. & Dreyfus G. (1993). Neural Networks and

Patan, K. (2008). *Artificial Neural Networks for the Modelling and Fault Diagnosis of Technical* 

Pati, Y. C. & Krishnaprasad, P. S. (1993). Analysis and synthesis of feedforward neural

Plett, G. L. (2003). Adaptive Inverse Control of Linear and Nonlinear Systems Using

Rivals, I. & Personnaz, L. (1996). Black Box Modeling With State-Space Neural Networks, in:

Roffel, B. & Betlem, B. (2006). *Process Dynamics and Control, Modeling for Control and* 

Saad Saoud, L. & Khellaf, A. (2009) Identification and Control of a Nonlinear Chemical

*Electrical Engineering Design and technologies*, Oct. 31- Nov. 2, 2009, Tunisia. Saad Saoud, L. & Khellaf, A. (2011) A Neural Network Based on an Inexpensive Eight Bit Microcontroller, *Neural computing and application*, 20(3), pp. 329-334. Saad Saoud, L.; Rahmoune, F.; Tourtchine, V. & Baddari, K. (2011). An Inexpensive

PhD Thesis, Yale University, New Haven, CT.

*Bulgarian Journal of Physics*. Vol. 31, pp. 163–169.

*Processes*, Springer-Verlag Berlin Heidelberg.

*Computation*, Vol.5 (2), pp. 165-199.

*Prediction*, John Wiley & Sons, 2006.

Accepted, July 26-28, 2011, Hangzhou , China.

*Networks,* vol. 4, pp. 73–85.

pp. 360-372.

Scientific, Singapore.

304.

pp.1119- 1132.

*ICAST* 2003.

USA, pp. 1614 – 1619.

asymptotically optimal approximators for functions of one variable, *in Proceeding IEEE International Conference on Neural Networks*, Orlando, FL, June 1994, pp. 299–

systems using neural networks," *International Journal of Control*, Vol. 71, No. 6,

neural network based cervical cancer diagnosis system: a microcontroller approach,

Predictive Controller, *American Control Conference,* 23-25 May 1990, San Diego, CA,

Nonlinear Adaptive Filtering: Unifying Concepts and New Algorithms, *Neural* 

networks using discrete affine wavelet transformations, *IEEE Transactions on Neural* 

Dynamic Neural Networks," IEEE *Transactions on Neural Networks*, Vol. 14, No. 2,

*Neural Adaptive Control Technology.* Zbikowski, R. & Hunt, K. J., pp. 237-264 World

process Plant Using Dynamical Neural Units, *Third International Conference on* 

Embedded Electronic Continuous Stirred Tank Reactor (CSTR) Based on Neural Networks, *The 2nd International Conference on Multimedia Technology (ICMT2011)*,

By a simple PIC and EEPROM programmer, we can charge the hexadecimal file (Wavelet.hex) and the binary file (CSTR.bin) which contain the hexadecimal file of the main assembly language program and data used for modeling and validation respectively. Verify the existence of the 9 V battery and that is all.

All necessary files can be found with the publisher of this book Intech.

#### **8. References**


By a simple PIC and EEPROM programmer, we can charge the hexadecimal file (Wavelet.hex) and the binary file (CSTR.bin) which contain the hexadecimal file of the main assembly language program and data used for modeling and validation respectively. Verify

Abiyev, R.H. (2003). Fuzzy Wavelet Neural Network for Control of Dynamic Plants,

Bakshi, B. R. & Stephanopoulos, G. (1993). Wave-net: A multi resolution hierarchical neural

Cotton, N. J.; Wilamowski, B. M. & Dündar, G. (2008). A Neural Network Implementation

Delyon, B.; Juditsky, A. & Benveniste, A. (1995). Accuracy analysis for wavelet approximations, *IEEE Transaction on Neural Networks*, vol. 6, pp. 332–348. De Moor, B. (1998). Daisy: Database for the identification of systems. Department of

http://www.esat.kuleuven.ac.be/sista/daisy/, Used data set: Continuous Stirred

Den Boer, W. (2005). *Active Matrix Liquid Crystal Displays: Fundamentals and Applications*.

Efe, M. Ö. (1996). Identification and Control of Nonlinear Dynamical Systems Using Neural

Efe, M. Ö. & Kaynak, O. (1997). Identification and Control of a Nonlinear Bioreactor Plant

Espinosa, J.; Vandewalle J. & Wertz V. (2005) *Fuzzy Logic, Identification and Predictive Control*,

Gao, X. (2002). A comparative research on wavelet neural networks. *Proceedings of the 9th* 

Gulbag, A.; Temurtas, F.; Tasaltin C. & Ozturk, Z. Z. (2009). A neural network implemented

Khalaf, N. M. A. (2006). Wavelet Network Identifier for Nonlinear Functions, thesis,

Henson, M. & Seborg D. (1990). Input-output linearization of general nonlinear processes,

Hong, J. (1992). Identification of stable systems by wavelet transform and artificial neural

Using Classical and Dynamical Neural Networks.' Proceeding of the International Symposium on Industrial Electronics (ISIE'97), Guimaraes, Portugal, Vol.3, 7-11

*International Conference on Neural Information Processing*, vol.4, pp. 1699 – 1703, 18-22

microcontroller system for quantitative classification of hazardous organic gases in the ambient air, *International Journal of Environment and Pollution*, Vol. 36-3, pp. 151 –

*Intelligent Engineering Systems*, February 25–29, 2008, Miami, Florida. Cybenko, G. (1989) Approximation by Superpositions of a Sigmoidal Function, *Mathematics* 

network with localized learning, *American Institute of Chemical Engineering Journal,*

on an Inexpensive Eight Bit Microcontroller, *12th International Conference on* 

*International Journal of Computational Intelligence*, 1 (2), pp.139-143.

the existence of the 9 V battery and that is all.

vol. 39, pp. 57–81.

Elsevier.

July 1997, pp. 1211-1215.

University of Technology, Iraq.

*AIChE Journal*, pp. 1753–1757.

networks, Ph.D. dissertation, Univ. Pittsburgh, PA.

Springer-Verlag, UK.

Nov. 2002.

165.

**8. References** 

All necessary files can be found with the publisher of this book Intech.

*of control, signals and systems*, Vol.2, pp. 303-314.

Networks. M.S. Thesis, Boðaziçi University.

Electrical Engineering -ESAT- K.U. Leuven, Belgium.

Tank Reactor, Section: Process Industry Systems, code: 98-002.


**Part 3** 

**Liquid Crystal Displays - Future Developments** 

