- Email: [email protected]
Comparison of two methods for oblique propagation in helicoidal bianisotropic mediums
Optics Communications 230 (2004) 369–386 www.elsevier.com/locate/optcom
Comparison of two methods for oblique propagation in helicoidal bianisotropic mediums John A. Polo Jr b
a,*
, Akhlesh Lakhtakia
b,1
a Department of Physics and Technology, Edinboro University of Pennsylvania, Edinboro, PA 16444, USA CATMAS – Computational & Theoretical Materials Sciences Group, Department of Engineering Science and Mechanics, Pennsylvania State University, University Park, PA 16802, USA
Received 13 September 2003; received in revised form 6 November 2003; accepted 12 November 2003
Abstract Two methods for oblique propagation in two dielectric examples of helicoidal bianisotropic mediums – chiral sculptured thin films and chiral smectic liquid crystals – are compared: (i) the series solution method, and (ii) the piecewise homogeneity approximation method. Both have to be numerically implemented, although only the former relies on an exact analytic solution of the frequency-domain Maxwell postulates. Convergence in the series solution method is investigated via calculations for a variety of planewave-illumination parameters and constitutive properties. Likewise, the maximum slice thickness for convergent results from the piecewise homogeneity approximation method is also studied. The computational efficiencies of both methods are compared. Ó 2003 Elsevier B.V. All rights reserved. PACS: 77.55.+f; 78.20.Ek; 78.20.Fm; 78.6.)w Keywords: Helicoidal bianisotropic medium; Matrix polynomial series; Propagation; Sculptured thin film; Smectic liquid crystal; Structural chirality
1. Introduction Chiral sculptured thin films (STFs) and chiral smectic liquid crystals (SLCs) represent two dielectric examples of helicoidal bianisotropic mediums (HBMs) in which the relative permittivity dyadic varies helicoidally along a fixed axis. Within a range of wavelengths – called the Bragg regime – that depends on the angle of incidence, circularly polarized light of one handedness is strongly reflected from these materials, but not circularly polarized light of the opposite handedness. The display of this discrimination is
*
1
Corresponding author. Tel.: +1-814-732-2655; fax: +1-814-732-2455. E-mail addresses: [email protected] (J.A. Polo Jr), [email protected] (A. Lakhtakia). Tel.: +1-814-863-4319; fax: +1-814-865-9974.
0030-4018/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2003.11.024
370
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
called the circular Bragg phenomenon. Parenthetically, another Bragg regime is also observable, at roughly twice the wavelength [1]; but as it does not display substantial discrimination between the two circular polarization states, it is technologically uninteresting. A host of literature exists on liquid crystals, as well as on the scientific and technological uses of the dependences of their optical responses on pressure, temperature and quasielectrostatic fields [2–4]. Chiral STFs can be tailor-made to a range of specifications to elicit particular optical responses. Consequently, there exists a large potential for practical applications of STFs, both as passive optical elements as well as optical sensors responding to mechanical and chemical environments; and some of that potential has already been actualized [1,5,6]. Exact closed-form solution of the frequency-domain Maxwell postulates in HBMs has been obtained only for propagation along the helicoidal axis (i.e., the axis of nonhomogeneity) [7–9], but not for oblique or off-axis propagation. An exact analytic solution for oblique propagation in the form of a matrix polynomial series became available in the last decade for cholesteric liquid crystals (CLCs) [10], which are the simplest of dielectric HBMs. The series converges uniformly although somewhat slowly, especially for large propagation distances. Elsewhere [11], 2 we have explored by numerical calculation the rate of convergence for CLCs over a range of constitutive parameters. The series solution method is also applicable to chiral STFs and chiral SLCs [1,12,13]. The present study is an extension of our initial exploration on the simple CLCs to these more complex HBMs. An algorithm to surmount the problem of convergence for large propagation distances is presented. The series solution method is compared with the piecewise homogeneity approximation method, when a plane wave is incident on a chiral STF/SLC of finite thickness. Although the latter method has been widely used [1,14–16], its computational efficiency remains unstudied. The plan of this paper is as follows. In Section 2, a 4 4 matrix ordinary differential equation for oblique propagation in a chiral STF/SLC is presented. The infinite series solution to this equation is recast in Section 3 in the form of a matrizant that has been related to spectral Green functions [17]. After presenting a technique to efficiently extend the series solution method to large propagation distances, the performance of the method for a model system is explored in Section 4. A useful portion of the parameter space encompassing the constitutive properties, the wavelength, and the angle of incidence is investigated. Section 5 contains a description of the piecewise homogeneity approximation method; and its performance on the previously chosen model system is investigated in Section 6. A comparison of the performances of both methods in Section 7 is followed by concluding remarks in Section 8.
2. Matrix ordinary differential equation The frequency-domain constitutive relations of a chiral STF/SLC can be written as [15] DðrÞ ¼ 0 S z ðzÞ S y ðvÞ oref S 1 S 1 ðzÞ EðrÞ ¼ ðzÞ EðrÞ; y z
ð1Þ
BðrÞ ¼ l0 H ðrÞ;
ð2Þ
where 0 and l0 are the permittivity and the permeability, respectively, of free space. As chiral STFs/SLCs are locally biaxial, the reference relative permittivity dyadic is given by oref ¼ a uz uz þ b ux ux þ c uy uy ;
2
~ 0 ðfÞ ¼ ½M ~ 0 ðnÞ½M ~ 0 ð1Þ‘ . Correction: Eq. (16) of this paper should read as follows: ½M
ð3Þ
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
371
where ux , uy and uz are the Cartesian unit vectors. The locally aciculate morphology of chiral STF/SLCs is described by the dyadic S y ðvÞ ¼ ðux ux þ uz uz Þ cos v þ ðuz ux ux uz Þ sin v þ uy uy ;
ð4Þ
where v 2 ½0; p=2 is the so-called angle of rise relative to the xy plane. The dyadic function S z ðzÞ ¼ ðux ux þ uy uy Þ cosðpz=XÞ þ ðuy ux ux uy Þ sinðpz=XÞ þ uz uz
ð5Þ
describes the helicoidal rotation of the permittivity dyadic ðzÞ, with 2X as the structural period along the z axis. The dielectric properties of a CLC are described as a special case by setting c ¼ a and v ¼ 0 [2,3]. In order to study the planewave response characteristics of a chiral STF/SLC, it is usual to confine the material between two infinite planes normal to the z axis. The wave vector of the incident plane wave lies at angle h with respect to the z axis, and at angle w with respect to the x axis in the xy plane. Provided the medium of incidence is vacuum, the quantity j ¼ k0 sin h is defined as a transverse wavenumber, with k0 being the wavenumber in vacuum. With the foregoing preliminaries, the constitutive relations (1) and (2) are substituted into the frequencydomain Maxwell curl postulates r EðrÞ ¼ ixl0 H ðrÞ;
r H ðrÞ ¼ ixðzÞ EðrÞ;
ð6Þ
where x is the angular frequency. A simplification in calculation results through the Oseen transformation [18] EðrÞ ¼ S z ðzÞ E0 ðrÞ;
H ðrÞ ¼ S z ðzÞ H 0 ðrÞ:
ð7Þ
Without loss of generality, the spatial Fourier transform of the fields is implemented with respect to variations in the xy plane; hence, the primed fields E0 ðrÞ and H 0 ðrÞ are expressed as E0 ðrÞ ¼ ½e0x ðz; j; wÞux þ e0y ðz; j; wÞuy þ e0z ðz; j; wÞuz exp½ijðx cos w þ y sin wÞ;
ð8Þ
H 0 ðrÞ ¼ ½h0x ðz; j; wÞux þ h0y ðz; j; wÞuy þ h0z ðz; j; wÞuz exp½ijðx cos w þ y sin wÞ;
ð9Þ
where w is an angle in the xy plane relative to the x axis, and j is a transverse wavenumber. Eq. (6) separate into two algebraic equations and four ordinary differential equations (ODEs), when (7)– (9) are substituted therein [12]. The two algebraic equations yield [1] j ðb a Þ 0 ½h0 ðz; j; wÞ sin n þ h0y ðz; j; wÞ cos n ex ðz; j; wÞ cos v sin v; x0 a s x a s j 0 h0z ðz; j; wÞ ¼ ½e ðz; j; wÞ sin n þ e0y ðz; j; wÞ cos n; xl0 x
e0z ðz; j; wÞ ¼
ð10Þ ð11Þ
where n¼
pz w; X
s ¼ cos2 v þ ðb =a Þ sin2 v:
ð12Þ
Thereafter, the four coupled ODEs can be represented compactly as the matrix ODE d 0 iX ½f ðf; j; wÞ ¼ ½P 0 ðf; j; wÞ ½f 0 ðf; j; wÞ; df p
ð13Þ
which has been written with the normalized axial distance f
pz X
ð14Þ
372
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
as the independent variable. The electromagnetic field components are contained in the column 4vector 2 0 3 ex ðf; j; wÞ 6 e0y ðf; j; wÞ 7 7 ð15Þ ½f 0 ðf; j; wÞ ¼ 6 4 h0x ðf; j; wÞ 5: 0 hy ðf; j; wÞ The 4 4 kernel matrix in (13) is 2 3 0 xl0 0 ip X 6 ip 0 xl0 0 7 X 7 ½P 0 ðf; j; wÞ ¼ 6 4 0 5 x0 c 0 ip X x0 b ip 0 0 X s 2 A cos n sin 2v 0 B sin n cos n 6 A sin n sin 2v 0 B sin2 n þ6 2 4 C sin n cos n C cos n 0 C sin2 n C sin n cos n A sin n sin 2v
3 B cos2 n B sin n cos n 7 7; 5 0 A cos n sin 2v
ð16Þ
where A¼
jðb a Þ ; 2a s
B¼
j2 ; x0 a s
C¼
j2 : xl0
ð17Þ
3. Series solution method Following Lakhtakia and Weiglhofer [12], the exact analytic (but not necessarily closed-form) solution of (13) can be written as the column-vector polynomial series 1 X ½f 0 ðf; j; wÞ ¼ ½pm ðj; wÞfm ; ð18Þ m¼0
which converges for all f because ½P 0 ðf; j; wÞ is a holomorphic function of f. For f ¼ 0, the polynomial series yields ½p0 ðj; wÞ ¼ ½f 0 ð0; j; wÞ; while the remaining coefficient 4-vectors on the right side of (18) are computed with the recurrence relation ½pm ðj; wÞ ¼
m1 1 iX X ½P 0 ðj; wÞ½pn ðj; wÞ; m p n¼0 mn1
m ¼ 1; 2; 3; . . .
The 4 4 matrixes in this recurrence relation are as follows: 2 3 0 xl0 0 ip X 6 ip 0 xl0 0 7 X 7 ½P 00 ðj; wÞ ¼ 6 4 0 5 x0 c 0 ip X x0 b ip 0 0 X s 2 3 B sin 2w B2 ½1 þ cos 2w A sin 2v cos w 0 2 B 6 A sin 2v sin w 0 ½1 cos 2w B2 sin 2w 7 2 7; þ6 C C 4 sin 2w 5 ½1 þ cos 2w 0 0 2 2 C C2 ½1 cos 2w sin 2w A sin 2v sin w A sin 2v cos w 2
ð19Þ
ð20Þ
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
3 A sin 2v sin w 0 B2m1 cos 2w B2m1 sin 2w imþ1 6 0 B2m1 sin 2w B2m1 cos 2w 7 6 A sin 2v cos w 7; ½P 0m ðj; wÞ ¼ m1 m1 4 5 0 0 m! C2 cos 2w C2 sin 2w m1 m1 C2 sin 2w C2 cos 2w A sin 2v cos w A sin 2v sin w
373
2
m ¼ 1; 3; 5; . . . ; ð21Þ
2
A sin 2v cos w m 6 i A sin 2v sin w ½P 0m ðj; wÞ ¼ 6 m! 4 C2m1 sin 2w C2m1 cos 2w
0 0 C2m1 cos 2w C2m1 sin 2w
B2m1 sin 2w B2m1 cos 2w 0 A sin 2v sin w
3 B2m1 cos 2w B2m1 sin 2w 7 7; 5 0 A sin 2v cos w
m ¼ 2; 4; 6; . . . ð22Þ 0
For solving boundary value problems, it is convenient to express the column 4-vector ½f ðf; j; wÞ in terms of a 4 4 matrix ½M 0 ðf; j; wÞ (called the matrizant), such that [1] ½f 0 ðf; j; wÞ ¼ ½M 0 ðf; j; wÞ½f 0 ð0; j; wÞ:
ð23Þ
It then follows that ½M 0 ðf; j; wÞ ¼
1 X
½R0m ðj; wÞfm ;
ð24Þ
m¼0
with ½R00 ðj; wÞ ¼ ½I; where ½I is the 4 4 identity matrix, and m1 1 iX X ½R0m ðj; wÞ ¼ ½P 0 ðj; wÞ½R0n ðj; wÞ; m p n¼0 mn1
ð25Þ
m ¼ 1; 2; 3; . . .
ð26Þ
4. Performance of series solution method A computer program, with double precision arithmetic used throughout, was written in FORTRAN 77 to calculate the matrizant ½M 0 ðf; j; wÞ. The program was executed on a desktop PC. In contrast to our earlier investigation on CLCs [11], in which terms were added in the series (24) for ½M 0 ðf; j; wÞ until all 16 elements reached a desired level of precision, the eight planewave remittances (i.e., four reflectances and four transmittances) of chiral STF/SLCs of finite thickness in the circularly polarized basis were calculated [1] and used to study convergence in the present study. It was found that, due to the complex nature of these dielectric HBMs in contrast to CLCs, the requirement that all elements in ½M 0 ðf; j; wÞ converge to a predetermined relative precision is unrealistic. In view of the use of these materials as optical elements, the planewave remittances serve as more practical measures of the accuracy of ½M 0 ðf; j; wÞ anyhow. Other target functions, such as the Jones matrix [2, Section 4.1.1], could have been used in lieu of the remittances; but we do not expect the general conclusions obtained here to be significantly altered. 4.1. Algorithm development for large f In the Bragg regime, the planewave reflectances saturate as the thickness of the chiral STF/SLC increases [15]. We chose the thickness as 60X – i.e., we computed the matrizant for f 60p – for the numerical
374
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
results reported here. This, however, presented a problem in calculating ½M 0 ðf; j; wÞ by the series solution method. The series (24) converges too slowly for f J p to be of any practical use for CLCs [11], but that problem is easily circumvented. Because the angle of rise is zero for CLCs, the f-period of the kernel matrix ½P 0 ðf; j; wÞ is p. Any f > 0 can be expressed as f ¼ ð‘ þ qÞp, where 0 < q < 1 and ‘ ¼ 0; 1; 2; . . . The matrizant for an arbitrary value of f can then be written as the product ½M 0 ðf; j; wÞ ¼ ½M 0 ðqp; j; wÞ½M 0 ðp; j; wÞ‘ . Thus, one never needs to directly use the series to calculate ½M 0 ðf; j; wÞ for f > p. The situation for chiral STFs/SLCs is not quite so simple. The f-period of the kernel matrix is 2p for j 6¼ 0, but only p when j ¼ 0 [19]. A compact algorithm should accommodate both conditions. Now, ½M 0 ðf; j; wÞ can be easily computed by the series (24) when 0 6 f 6 p, but not when p < f 6 2p. In other words, if we rewrite ½M 0 ðp þ Df; j; wÞ ¼ ½Q0 ðDf; j; wÞ½M 0 ðp; j; wÞ;
0 < Df 6 p;
ð27Þ
the key issue is the computation of the auxiliary matrizant ½Q0 ðDf; j; wÞ. An expression for the auxiliary matrizant can be derived as follows: After a displacement of f ¼ p, the permittivity dyadic ðfÞ is rotated 180° about the z axis. If the coordinate system is also rotated by 180° about the z axis, the elements of ðfÞ take on the same numerical values in the second half-period as in the first half-period. In the rotated coordinate system however, w becomes p þ w. Therefore, ½Q0 ðDf; j; wÞ ¼ ½I ½M 0 ðDf; j; p þ wÞ ½I ¼ ½M 0 ðDf; j; p þ wÞ; 0 < Df 6 p; ð28Þ where the two factors of ð½IÞ represent a similarity transformation. Any f > 0 can be represented as f ¼ ð2j þ ‘ þ qÞp, where 0 < q < 1; either ‘ ¼ 0 or ‘ ¼ 1; and j ¼ 0; 1; 2; . . . Therefore, 8 j > < ½M 0 ðqp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ; ‘ ¼ 0; j ð29Þ ½M 0 ðf; j; wÞ ¼ > : ½M 0 ðqp; j; p þ wÞ½M 0 ðp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ; ‘ ¼ 1: The two conditional parts of the foregoing equation can be combined to yield the equation j ‘ ½M 0 ðf; j; wÞ ¼ ½M 0 ðqp; j; ‘p þ wÞ½M 0 ðp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ;
ð30Þ
which holds for both ‘ ¼ 0 and ‘ ¼ 1. Thus, the matrizant can be computed using the series (24) for any value of the normalized axial distance f. 4.2. Model system The number of terms in the series (24) that can be summed is limited to 170 on our desktop PC, due to overflow in calculating m!. Therefore, the 170-term sum was used as a best estimate of the true matrizant by which the convergence of the series (24) could be judged. To characterize the performance of the series solution method, the summation was truncated at m ¼ N , and the value Nreq of N required to yield all remittances to within 0.1% of the values predicted by the 170-term calculation of ½M 0 ðf; j; wÞ was determined for a range of parameters describing both the incident plane wave and the constitutive properties. The starting model system for the present study is delineated by: X ¼ 200 nm, a ¼ 2, b ¼ 3, c ¼ 2:4, and v ¼ 25°. The wavelength of the incident light in free space is k0 ¼ 600 nm and the plane of incidence is described by w ¼ 0°. Except for the determination of Nreq vs. f, the thickness of the chiral STF/SLC was set
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
375
equal to 61:5X. This value was chosen so that not only was there saturation in the values of the remittances with respect to f [1], but also that our algorithm was tested on all three classes of matrixes appearing on the right side of (30). The transverse wavenumber j ¼ ð2p=k0 Þ sin h, where h is the angle of planewave incidence with respect to the z axis. Axial propagation occurs inside the chiral STF/SLC for h ¼ 0°, while oblique propagation requires h 2 ð0°; 90°Þ. 4.3. Nreq vs. f Fig. 1 shows the curves of Nreq vs. f for h ¼ 0°; 30°; 60°, 89°, and 0 < f K 200. Calculation of the set of matrizant factors in (30) has an algorithmic periodicity (not to be confused with the structural periodicity) in f of 4p. However, as f increases, the last factor in (30) is raised to ever higher powers j – thereby leading to the conjecture that, in order to maintain the same level of accuracy in the final result, higher precision in this factor might be required as f is increased. This does not appear to be the case. After a very sharp rise in Nreq for small f, the curves almost level off. A small amount of oscillation is observed in every curve. It is most pronounced at h ¼ 60° and amounts to approximately a 10% fluctuation. The maximum value of Nreq for all h chosen is slightly greater than 65. Parenthetically, the noisiness of data in Fig. 1 (as also Figs. 2–5 and 8–11) can be attributed to two reasons: (i) the periodic nonhomogeneity of HBMs, which is responsible for the differential transmissivity characteristic of the circular Bragg phenomenon; and (ii) because the quantities plotted either vary or were allowed to vary in discrete steps. 4.4. Nreq vs. X Fig. 2 shows the curves of Nreq vs. X 2 ½50; 600 nm, for h ¼ 0°; 30°; 60°, and 89°. The curves let us conclude that Nreq X, which can be explained in the following manner: As the period 2X increases, the morphological length-scale becomes comparable to the wavelength, thereby enhancing the complexity of the spatial distribution of the electromagnetic field. The increased complexity entails a larger number of
70
60
Nreq
50
40
30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10
0 0
50
100
150
200
ζ Fig. 1. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the normalized axial distance f ¼ pz=X for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and v ¼ 25°.
376
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 160 θ
140
0 deg. 30 deg. 60 deg. 89 deg.
120
Nreq
100 80 60 40 20 0 0
200
400
600
Fig. 2. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the half-period X for: k0 ¼ 600 nm, w ¼ 0°, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
terms in the series (24) for convergence. Additionally, we note that Nreq for h ¼ 89° is a little more than three times that for h ¼ 0°. 4.5. Nreq vs. c Fig. 3 displays the curves of Nreq vs. c , with c ranging from considerably below a to considerably above a . These curves do not evince an easily identifiable pattern. The curves for lower values of h are
100
θ 0 deg. 30 deg. 60 deg. 89 deg.
90 80 70
Nreq
60 50 40 30 20 10 0 1.5
2
εc
2.5
3
Fig. 3. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus c for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
377
nearly flat. Those for higher values of h show a minimum in Nreq at c ¼ 2:7, and the depth of this minimum increases with the angle of incidence. As with the other sets of calculations reported here, Nreq rises with increasing h. 4.6. Nreq vs. v The relationship between Nreq and v, shown in Fig. 4, is nearly flat. There is one notable exception at v ¼ 45° where there is a very sharp peak in the h ¼ 0° curve. This peak corresponds to the pseudoisotropic point (i.e., c ¼ b =s) [20]. At the pseudoisotropic point, the Bragg regime vanishes identically for normal incidence and is seriously weakened for oblique incidence. The general trend of larger Nreq for larger h is also evident in Fig. 4. 4.7. Nreq vs. precision The value of Nreq was calculated for relative error in the remittances ranging four orders of magnitude from 105 to 101 . The results of these calculations are presented in Fig. 5. For the chosen values of h, the maximum change in Nreq over the large range of relative precision in all remittances is about 60%. 4.8. Calculation time The calculations presented here were carried out on a desktop PC with a clock speed of 1.8 GHz. The time required to calculate the matrizant was determined by timing 100 iterations of the calculation. This number of iterations gave consistent results when the calculation time Tseries (in millisecond) was remeasured, but was not so long as to be tiresome. The time to complete the calculation depends only on the number of terms calculated in the series. The variation of Tseries with respect to Nreq is displayed in Fig. 6. The data were fit with a quadratic with excellent results; indeed, on the scale of the figure, the collected data and the quadratic fit
80 70 60
Nreq
50 40 30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10 0 0
50
χ (deg.)
100
Fig. 4. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the rise angle v for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and z ¼ 61:5X.
378
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 80 70 60
Nreq
50 40 30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10 0
10-5
10-4
10-3 Precision
10-2
1
10-1
Fig. 5. Series solution method. Number of terms required for various values of relative precision in remittances for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
250
Tseries (ms)
200
150
100
50
0 0
50
100
150
200
Nreq Fig. 6. Series solution method. Calculation time versus the number of terms.
2 Tseries ¼ 8 103 Nreq þ 7:6 102 Nreq
ð31Þ
are indistinguishable. Considering the way in which the matrix coefficients of the series (24) are calculated recursively as per (26), we might have anticipated the quadratic time-dependence of the series solution method. For the coefficient matrix of the mth term in the series, m matrix multiplications are required on the right side of (26). A total of ð1=2ÞNreq ðNreq þ 1Þ matrix multiplications must be carried out for a series with 2 Nreq terms. For large Nreq , the number of matrix multiplications is thus proportional to Nreq . Except for very large values of X, the number of terms required for all calculations reported here was Nreq < 90. Therefore, the calculation time by the series solution method was Tseries 6 72 ms.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
379
5. Piecewise homogeneity approximation method The piecewise homogeneity approximation method has been successfully applied to chiral STFs/SLCs [1,14]. In this method, the material – which is nonhomogeneous along the z axis – is divided into a series of slices (or pieces) perpendicular to that axis, each of thickness Dzsl . Each slice is treated as if it were homogenous with its permittivity dyadic taken to be that at the center of the slice. In a medium with uniform permittivity dyadic, ½P 0 ðf; j; XÞ becomes independent of f [9], and the solution to (13) can be written in closed form as [21, Section 2.5] iX 0 ½P ðj; wÞDf ½f ðf; j; wÞ; Df ¼ f0 f: ½f 0 ðf0 ; j; wÞ ¼ exp ð32Þ p Therefore, if fj1 represents the normalized axial distance of the start of the jth slice, and Dfsl is the normalized slice thickness, the transfer of electromagnetic fields across this slice is expressed by iX 0 Df ½P fj1 þ sl ; j; w Dfsl ½f ðfj1 ; j; wÞ; j ¼ 1; 2; 3; . . . ½f 0 ðfj ; j; wÞ ¼ exp ð33Þ p 2 Denoting the matrizant describing propagation through the jth slice by iX 0 Df ½M 0j ðj; wÞ ¼ exp ½P fj1 þ sl ; j; w Dfsl ; p 2
ð34Þ
we see that propagation through a normalized axial distance f ¼ Nsl Dfsl is described by ½M 0 ðf; j; wÞ ½M 0N ðj; wÞ½M 0N sl
sl 1
ðj; wÞ ½M 01 ðj; wÞ;
ð35Þ
where the number of slices Nsl must be adequate for a convergent solution. In order to execute the calculation of ½M 0 ðf; j; wÞ via (35), the matrix exponential in (34) must be determined. Now, as n 1 X iX 0 iX n ½P ðfj1 þ Dfsl =2; j; wÞDfsl ¼ ½I þ Dfsl ½P 0 ðfj1 þ Dfsl =2; j; wÞ exp ð36Þ p p n¼1 is an infinite series, the sum matrix on the right side of (36) was computed term by term. Following Venugopal and Lakhtakia [1], the summation was stopped at that value of n for which the changes in all elements of the sum matrix were less than or equal to one part in 1012 . Parenthetically, we note that the left side of (36) can be evaluated via diagonalization of ½P 0 ðfj1 þ Dfsl =2; j; wÞ, following eigenvalue analysis [22], provided no eigenvalue has a geometric multiplicity other than unity [21, Section 2.5].
6. Performance of piecewise homogeneity approximation method The accuracy of the piecewise homogeneity approximation method is controlled by Dfsl . However, there is no general way of knowing, in advance, how thin to make the slices in order to obtain a particular accuracy in the calculation. To obtain a desired degree of accuracy, one must repeat the calculation several times, increasing the number of slices each time, and watch for stabilization of the remittances. We carried out a series of calculations to look at how the required Dfsl varies as a function of the parameters describing the incident plane wave and the chiral STF/SLC. The accuracy of the approximation was judged by comparing results to those of the 170-term series solution of Section 4: all remittances had to be within 0.1% of those from the series solution. The parameters and the range of independent variables studied are the same as those for studying the series solution method.
380
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
6.1. Algorithm development Although there is no difficulty in calculating the entire matrizant directly using the piecewise homogeneity approximation method for any f, the calculation time for large f can be reduced substantially by exploiting the periodicity of chiral STFs/SLCs. The factorization is simpler than for the series solution method, since the piecewise computation evinces no problems for ½M 0 ð2p; j; wÞ. Any f > 0 can be expressed as f ¼ 2ð‘ þ qÞp, where 0 < q < 1 and ‘ ¼ 0; 1; 2; . . . The matrizant for an arbitrary value of f can then be obtained as the product ‘
½M 0 ðf; j; wÞ ¼ ½M 0 ð2qp; j; wÞ½M 0 ð2p; j; wÞ :
ð37Þ
In order to determine Dfsl , Nsl was increased until the desired accuracy in the remittances was obtained. The number of slices was partitioned between the two intervals 0 6 f 6 2qp and 0 6 f 6 2p – associated with the two factors on the right side of (37) – in proportion to the size of each interval and rounded to the nearest integer. The average slice thickness was then taken to be hDzsl i ¼ ð1 þ qÞ
2X ; Nsl
ð38Þ
which quantityÕs variation with various incidence and constitutive parameters was examined. 6.2. hDzsl i vs. h Fig. 7 shows hDzsl i vs. h for the starting model system of Section 4.2. Most notable is the very large spike at h ¼ 0°. This is to be expected since it is possible to calculate ½M 0 ðf; j; wÞ for h ¼ 0° with just one slice for any f [15]. As we chose f ¼ 61:5p and our algorithm is based on (37), at least two slices are required; accordingly, hDzsl i ¼ 350 nm for X ¼ 200 nm, just as found. Because hDzsl i for h ¼ 0° is about two orders of magnitude greater than for h only slightly greater than 0° for all studies reported here, the curves for h ¼ 0° are not shown in the subsequent figures in order to not obscure the results for other values of h.
400 350
<∆Zsl> (nm)
300 250 200 150 100 50 0 0
10
θ (deg.)
20
Fig. 7. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus angle of incidence h for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
381
6.3. hDzsl i vs. f Fig. 8(a) and (b) shows the variation of hDzsl i with f for h ¼ 30°, 60° and 89°. The h ¼ 30° curve in Fig. 8(a) is essentially flat. The h ¼ 89° curve shows an initial steep decline tapering off to a much more gradual decline as a function f, while the curve for h ¼ 60° in Fig. 8(b) shows a nearly steady decline. This greater effort required to calculate the matrizant (smaller slice thickness) as f increases, is contrary to what was observed for the series solution with Nreq relatively constant. Furthermore, all three curves have a significant amount of fluctuation about the general trend, with the h ¼ 60° curve showing the largest amount. hDzsl i changes by approximately a factor of 5 in one cycle in the most extreme fluctuation. This is more than an order of magnitude greater fluctuation than the 10% effect seen in Nreq in the series solution. Over the entire range of calculations, hDzsl i changes by greater than a factor of 20. In contrast, the change in Nreq for the series solution, excluding a few points at f < 1, changes only by roughly a factor of 3 over the entire range of calculation. These observations allow us to conclude that the error in replacing a thin nonhomogeneous slice by a homogeneous slice in the piecewise homogeneity approximation method accumulates as the normalized axial thickness f increases. In order to offset the accumulation of error, hDzsl i has to be reduced. 6.4. hDzsl i vs. X The variation of hDzsl i with X is shown in Fig. 9 for h ¼ 30°, 60°, 89°. For h ¼ 60° and 89°, there is a general tendency to larger hDzsl i as X increases. The h ¼ 30° curve shows a similar trend but with a broad dip near X ¼ 500 nm, the dip corresponding to virtually null values of the cross-polarized transmittances and the co-polarized reflectances. Leaving aside that broad dip, we conclude that the number of slices increases with X. This is in line with the series solution method for which Nreq (and hence the computational effort) increases with X for all angles of incidence. For both methods, the computational effort is more at higher h.
7
7 θ
6
6 θ
30 deg.
4 3
4 3
2
2
1
1 0
0 0
(a)
60 deg.
5
89 deg.
<∆Zsl> (nm)
<∆Zsl> (nm)
5
50
100
ζ
150
0
200
(b)
50
100
ζ
150
200
Fig. 8. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the normalized axial distance f ¼ pz=X for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°.
382
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 10 9
θ
8
30 deg. 60 deg. 89 deg.
<∆Zsl> (nm)
7 6 5 4 3 2 1 0 0
200
400
600
Ω (nm)
Fig. 9. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the half-period X for: k0 ¼ 600 nm, w ¼ 0°, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
6.5. hDzsl i vs. c There is no general trend evinced by the hDzsl i vs. c curves in Fig. 10. Oscillations in these curves, however, are quite dramatic – most prominently for h ¼ 30° when c K 2:25. The oscillations are in marked contrast to the dependence of Nreq on c in Fig. 3 for the series solution method. The oscillations in the hDzsl i vs. c curve coincide with oscillations in the remittances inside the Bragg regime. Furthermore, the drop for h ¼ 30° when c J 2:25 corresponds to drops in i(i) the co-polarized reflectance when the incident plane wave is right circularly polarized, and (ii) the co-polarized transmittance when the incident plane wave is left circularly polarized, 9 θ
8
30 deg. 60 deg. 89 deg.
7
<∆Zsl> (nm)
6 5 4 3 2 1 0 1.5
2.0
εc
2.5
3.0
Fig. 10. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus c for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
383
outside the Bragg regime. The correspondence between the oscillations in Fig. 10 and the remittances is less striking when h ¼ 60°, and virtually vanishes when h ¼ 89°. The hDzsl i-range for calculations at all h is quite large. The smallest value of hDzsl i in the h ¼ 89° curve is 0.0772 nm, while the largest value in the h ¼ 30° curve is 8.14 nm. The two values differ by slightly more than two order of magnitude. This is very much greater than the range of values seen in Nreq of the series solution for which the maximum value of Nreq differs from the minimum only by a factor of roughly 4. 6.6. hDzsl i vs. v The curves of hDzsl i vs. v are shown in Fig. 11 for h ¼ 30°, 60°, and 89°. Again, the fluctuations in these curves are extremely pronounced. Some, but not all, features of the hDzsl i vs. v curves can be correlated to remittance variations. Recall that, for the series method, Nreq vs. v has only one peak at v ¼ 45° (the pseudoisotropic point) when h ¼ 0°. Over the range of calculations illustrated in Figs. 4 and 11, hDzsl i changes by a factor of nearly 47, while the maximum value of Nreq differs from the minimum by less than a factor of 4. 6.7. hDzsl i vs. precision Fig. 12(a) shows the relationship between hDzsl i and relative precision in the remittances which range from 101 to 105 . Unlike the series solution method – for which Nreq varies very little with relative precision (see Fig. 5) – hDzsl i depends very strongly on relative precision. The log–log plot in Fig. 12(b) shows that hDzsl i varies as the square of the relative precision. From a relative precision of 105 to 101 , hDzsl i changes by approximately two orders of magnitude; in contrast, the maximum change in Nreq is by 60% in Fig. 5. 6.8. Calculation time Calculation times with the piecewise homogeneity approximation method were determined in the same way as for the series solution method. There are two main factors which influence the calculation time Tpiece 8 θ
7
30 deg. 60 deg. 89 deg.
<∆zsl> (nm)
6 5 4 3 2 1 0 0
20
40
χ (deg.)
608
0
Fig. 11. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the rise angle v for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and z ¼ 61:5X.
384
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
30 10 1
20
<∆Zsl> (nm)
<∆Zsl> (nm)
θ 30 deg. 60 deg. 89 deg.
10 0 θ
10
0 -5 10
(a)
30 deg. 60 deg. 89 deg.
10-1
10- 4
10- 3
10- 2
10-2
10- 1
10-5
10-4
(b)
Precision
10-3
10-2
10-1
Precision
Fig. 12. Piecewise homogeneity approximation method. Slice thickness required for various values of relative precision in remittances for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
for the piecewise method: X and hDzsl i. Measurements of Tpiece were made for hDzsl i ¼ 0:1, 0.5, 1.0, 5.0, and 10.0 nm, with X ranging from 50 nm to 600 nm. The results for hDzsl i ¼ 0:1, 1.0, and 5.0 nm are reported in Fig. 13. As can be deduced from this figure, Tpiece varies linearly with X for all values of hDzsl i. The dependences of the slopes of the Tpiece vs. X curves on hDzsl i can be represented by a simple power law. A reasonable approximation for calculation time emerges from Fig. 13 as Tpiece ¼ 0:20XhDzsl i0:89 ;
ð39Þ
where Tpiece is in millisecond and hDzsl i in nm. This formula predicts the calculation time within 10% for most cases, although the error may be as much as 25% at low values of X and high values of hDzsl i.
1000 <∆ Zsl> 0.1 nm 1.0 nm 5.0 nm
800
Tpiece(ms)
600
400
200
0 0
200
400
Ω (nm)
600
Fig. 13. Piecewise homogeneity approximation method. Calculation time versus X for various values of hDzsl i.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
385
7. Performance comparison The series solution method and the piecewise uniform permittivity method offer two viable means of determining the propagation of light in chiral STFs and chiral SLCs. The former method is exact and analytic, whereas the latter method is heuristic. For k0 > 3X and the chosen values of the constitutive parameters, planewave remittances can be calculated to a precision of 0:1% with less than 90 terms in the series (24). With this number of terms, calculations can easily be performed on a desktop computer, and a single calculation of ½M 0 ðf; j; wÞ appears instantaneous to the user with a calculation time 672 ms. The number of terms required for higher precision is only modestly higher, still well within the range accessible by desktop computers. To increase the relative precision from 0:1% to 0:001%, the number of terms in (24) only needs to be increased by approximately 30% in the worst case that we studied, resulting in a 70% increase in calculation time. The required number of terms Nreq for a given accuracy is very weakly dependant on the parameters describing the material and the incident plane wave. A general trend is that Nreq increases with angle of incidence h. These conclusions are true both inside as well as outside the Bragg regime. With the exception of the h ¼ 89°, all Nreq vs. f curves in Fig. 1 pass through the Bragg regime. There is no substantial change evident in Nreq as this occurs. In Figs. 2–5, the h ¼ 0° and h ¼ 30° curves are in the Bragg regime, while the curves for other values of h are not. No qualitative difference is apparent between those curves that are in the Bragg regime, and those that are not. The performance of the piecewise homogeneous permittivity approximation method appears to be much more variable. The slice thickness required to obtain 0.1% precision in remittances fluctuates dramatically. For all calculations made with 0.1% precision and X ¼ 200 nm, hDzsl i varies between 0.0772 nm and 8.14 nm. This implies that the calculation time varies roughly between 6 ms and 400 ms. To increase relative precision form 0.1% to 0.001%, hDzsl i must be reduced about an order of magnitude, increasing the calculation time by nearly a factor of 8. The calculation times with both the series solution method and the piecewise homogeneity approximation method vary with the parameters describing the helicoidal bianisotropic medium and the incident plane wave. A general conclusion is that calculation times for the former method vary less erratically than for the latter method. For the following discussion, all calculation times were obtained using (31) and (39). The calculation time with either method increases with the angle of incidence of the plane wave. With the piecewise homogeneity approximation method, there is a large change in the first few degrees away from normal incidence, amounting to nearly an order of magnitude. For our starting model system, the piecewise homogeneity approximation method is roughly 10 times faster than the series solution method at h ¼ 0°; the two methods are comparable for h J 2°; but the piecewise homogeneity approximation method performs at about half the speed of the series solution method at h ¼ 89°. The calculation time with either method increases with the structural period. The rate of increase of the calculation time with the structural period is larger for the series solution method than for the piecewise approximation method. For the model system, the calculation times for the two methods are very similar when X ¼ 50 nm, but the piecewise homogeneity approximation method is faster (except at h ¼ 89°) when X ¼ 600 nm. Thus, for h ¼ 60°, Tpiece ¼ 12 ms at X ¼ 50 nm and 45 ms at X ¼ 600 nm; whereas Tseries ¼ 9 and 120 ms, respectively. Over the range of X studied, neither method appears to have a clear advantage over the other. The calculation times for the two methods clearly differ as functions of f. Except at extremely small f, Tseries is independent of f. This is true of the piecewise homogeneity approximation method at h ¼ 30° as well, but at larger angles Tpiece increases with f. At h ¼ 89°, for instance, from f ¼ 5 to f ¼ 190, Tpiece changes from about 20 to 100 ms. The series solution method, on the other hand, has a roughly constant calculation time of about 35 ms.
386
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
No general trends regarding the calculation times emerged with respect to either the ratio c =a or v, except that the larger fluctuations in hDzsl i than in Nreq imply that Tpiece varies more erratically than Tseries with respect to these parameters.
8. Concluding remarks The exact and analytic series solution method can be applied only if the 4 4 kernel matrix in (13) is a holomorphic function of f, while the heuristic piecewise homogeneity approximation method is definitely more versatile. However, as we have shown here, the former method does offer some advantages when determining the oblique propagation characteristics of light in chiral STFs and chiral SLCs. The number of terms required in the series (24) to calculate planewave remittances to a given precision remains relatively constant over a wide range of parameters. The same cannot be said of the average slice thickness needed to implement the piecewise homogeneity approximation method. Also, to calculate the remittances for a specific set of geometric and constitutive parameters to a predetermined level of accuracy, the entire piecewise calculation has to be repeated several times with different values of Nsl . Thus, the actual calculation time for a single calculation with the piecewise homogeneity approximation method could be significantly larger than those reported here for which the slice thicknesses were pre-determined. In contrast, the convergence of the series solution method can be monitored as the series is summed and the process terminated when the desired accuracy has been achieved. This could be a significant advantage when a large set of calculations is to be performed.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
V.C. Venugopal, A. Lakhtakia, Proc. R. Soc. Lond. A 456 (2000) 125. S. Chandrasekhar, Liquid Crystals, second ed., Cambridge University Press, Cambridge, UK, 1992. P.G. de Gennes, J. Prost, The Physics of Liquid Crystals, second ed., Clarendon Press, Oxford, UK, 1993. S.D. Jacobs (Ed.), Selected Papers on Liquid Crystals for Optics, SPIE Press, Bellingham, WA, 1992. I.J. Hodgkinson, Q.h. Wu, Adv. Mater. 13 (2001) 889. A. Lakhtakia, Mater. Sci. Eng. C 19 (2002) 427. E.I. Kats, Sov. Phys. JETP 32 (1971) 1004. R. Nityananda, Mol. Cryst. Liq. Cryst. 21 (1973) 315. A. Lakhtakia, W.S. Weiglhofer, Microw. Opt. Technol. Lett. 6 (1993) 804. A. Lakhtakia, W.S. Weiglhofer, Microw. Opt. Technol. Lett. 12 (1996) 245. J.A. Polo Jr., A. Lakhtakia, Microw. Opt. Technol. Lett. 35 (2002) 397. A. Lakhtakia, W.S. Weiglhofer, Proc. R. Soc. Lond. A 453 (1997) 93, corrections: 454 (1998) 3275. M. Schubert, C.M. Herzinger, Phys. Stat. Sol. (a) 188 (2001) 1563. I. Abdulhalim, R. Weil, L. Benguigui, Liq. Cryst. 1 (1986) 155. V.C. Venugopal, A. Lakhtakia, in: O.N. Singh, A. Lakhtakia (Eds.), Electromagnetic Fields in Unconventional Materials and Structures, Wiley, New York, NY, 2000 (Chapter 5). J.A. Sherwin, A. Lakhtakia, Opt. Commun. 222 (2003) 305. A. Lakhtakia, W.S. Weiglhofer, IEE Proc.-Microw. Antennas Propag. 144 (1997) 57. C.W. Oseen, J. Chem. Soc. Faraday Trans. II 29 (1933) 883. € bertrag. 53 (1999) 287. A. Lakhtakia, V.C. Venugopal, Arch. Elektron. U A. Lakhtakia, Microw. Opt. Technol. Lett. 34 (2002) 367. H. Hochstadt, Differential Equations: A Modern Approach, Dover Press, New York, NY, 1975 (Section 2.5). M. Schubert, Phys. Rev. B 53 (1996) 4265.
Comparison of two methods for oblique propagation in helicoidal bianisotropic mediums John A. Polo Jr b
a,*
, Akhlesh Lakhtakia
b,1
a Department of Physics and Technology, Edinboro University of Pennsylvania, Edinboro, PA 16444, USA CATMAS – Computational & Theoretical Materials Sciences Group, Department of Engineering Science and Mechanics, Pennsylvania State University, University Park, PA 16802, USA
Received 13 September 2003; received in revised form 6 November 2003; accepted 12 November 2003
Abstract Two methods for oblique propagation in two dielectric examples of helicoidal bianisotropic mediums – chiral sculptured thin films and chiral smectic liquid crystals – are compared: (i) the series solution method, and (ii) the piecewise homogeneity approximation method. Both have to be numerically implemented, although only the former relies on an exact analytic solution of the frequency-domain Maxwell postulates. Convergence in the series solution method is investigated via calculations for a variety of planewave-illumination parameters and constitutive properties. Likewise, the maximum slice thickness for convergent results from the piecewise homogeneity approximation method is also studied. The computational efficiencies of both methods are compared. Ó 2003 Elsevier B.V. All rights reserved. PACS: 77.55.+f; 78.20.Ek; 78.20.Fm; 78.6.)w Keywords: Helicoidal bianisotropic medium; Matrix polynomial series; Propagation; Sculptured thin film; Smectic liquid crystal; Structural chirality
1. Introduction Chiral sculptured thin films (STFs) and chiral smectic liquid crystals (SLCs) represent two dielectric examples of helicoidal bianisotropic mediums (HBMs) in which the relative permittivity dyadic varies helicoidally along a fixed axis. Within a range of wavelengths – called the Bragg regime – that depends on the angle of incidence, circularly polarized light of one handedness is strongly reflected from these materials, but not circularly polarized light of the opposite handedness. The display of this discrimination is
*
1
Corresponding author. Tel.: +1-814-732-2655; fax: +1-814-732-2455. E-mail addresses: [email protected] (J.A. Polo Jr), [email protected] (A. Lakhtakia). Tel.: +1-814-863-4319; fax: +1-814-865-9974.
0030-4018/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2003.11.024
370
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
called the circular Bragg phenomenon. Parenthetically, another Bragg regime is also observable, at roughly twice the wavelength [1]; but as it does not display substantial discrimination between the two circular polarization states, it is technologically uninteresting. A host of literature exists on liquid crystals, as well as on the scientific and technological uses of the dependences of their optical responses on pressure, temperature and quasielectrostatic fields [2–4]. Chiral STFs can be tailor-made to a range of specifications to elicit particular optical responses. Consequently, there exists a large potential for practical applications of STFs, both as passive optical elements as well as optical sensors responding to mechanical and chemical environments; and some of that potential has already been actualized [1,5,6]. Exact closed-form solution of the frequency-domain Maxwell postulates in HBMs has been obtained only for propagation along the helicoidal axis (i.e., the axis of nonhomogeneity) [7–9], but not for oblique or off-axis propagation. An exact analytic solution for oblique propagation in the form of a matrix polynomial series became available in the last decade for cholesteric liquid crystals (CLCs) [10], which are the simplest of dielectric HBMs. The series converges uniformly although somewhat slowly, especially for large propagation distances. Elsewhere [11], 2 we have explored by numerical calculation the rate of convergence for CLCs over a range of constitutive parameters. The series solution method is also applicable to chiral STFs and chiral SLCs [1,12,13]. The present study is an extension of our initial exploration on the simple CLCs to these more complex HBMs. An algorithm to surmount the problem of convergence for large propagation distances is presented. The series solution method is compared with the piecewise homogeneity approximation method, when a plane wave is incident on a chiral STF/SLC of finite thickness. Although the latter method has been widely used [1,14–16], its computational efficiency remains unstudied. The plan of this paper is as follows. In Section 2, a 4 4 matrix ordinary differential equation for oblique propagation in a chiral STF/SLC is presented. The infinite series solution to this equation is recast in Section 3 in the form of a matrizant that has been related to spectral Green functions [17]. After presenting a technique to efficiently extend the series solution method to large propagation distances, the performance of the method for a model system is explored in Section 4. A useful portion of the parameter space encompassing the constitutive properties, the wavelength, and the angle of incidence is investigated. Section 5 contains a description of the piecewise homogeneity approximation method; and its performance on the previously chosen model system is investigated in Section 6. A comparison of the performances of both methods in Section 7 is followed by concluding remarks in Section 8.
2. Matrix ordinary differential equation The frequency-domain constitutive relations of a chiral STF/SLC can be written as [15] DðrÞ ¼ 0 S z ðzÞ S y ðvÞ oref S 1 S 1 ðzÞ EðrÞ ¼ ðzÞ EðrÞ; y z
ð1Þ
BðrÞ ¼ l0 H ðrÞ;
ð2Þ
where 0 and l0 are the permittivity and the permeability, respectively, of free space. As chiral STFs/SLCs are locally biaxial, the reference relative permittivity dyadic is given by oref ¼ a uz uz þ b ux ux þ c uy uy ;
2
~ 0 ðfÞ ¼ ½M ~ 0 ðnÞ½M ~ 0 ð1Þ‘ . Correction: Eq. (16) of this paper should read as follows: ½M
ð3Þ
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
371
where ux , uy and uz are the Cartesian unit vectors. The locally aciculate morphology of chiral STF/SLCs is described by the dyadic S y ðvÞ ¼ ðux ux þ uz uz Þ cos v þ ðuz ux ux uz Þ sin v þ uy uy ;
ð4Þ
where v 2 ½0; p=2 is the so-called angle of rise relative to the xy plane. The dyadic function S z ðzÞ ¼ ðux ux þ uy uy Þ cosðpz=XÞ þ ðuy ux ux uy Þ sinðpz=XÞ þ uz uz
ð5Þ
describes the helicoidal rotation of the permittivity dyadic ðzÞ, with 2X as the structural period along the z axis. The dielectric properties of a CLC are described as a special case by setting c ¼ a and v ¼ 0 [2,3]. In order to study the planewave response characteristics of a chiral STF/SLC, it is usual to confine the material between two infinite planes normal to the z axis. The wave vector of the incident plane wave lies at angle h with respect to the z axis, and at angle w with respect to the x axis in the xy plane. Provided the medium of incidence is vacuum, the quantity j ¼ k0 sin h is defined as a transverse wavenumber, with k0 being the wavenumber in vacuum. With the foregoing preliminaries, the constitutive relations (1) and (2) are substituted into the frequencydomain Maxwell curl postulates r EðrÞ ¼ ixl0 H ðrÞ;
r H ðrÞ ¼ ixðzÞ EðrÞ;
ð6Þ
where x is the angular frequency. A simplification in calculation results through the Oseen transformation [18] EðrÞ ¼ S z ðzÞ E0 ðrÞ;
H ðrÞ ¼ S z ðzÞ H 0 ðrÞ:
ð7Þ
Without loss of generality, the spatial Fourier transform of the fields is implemented with respect to variations in the xy plane; hence, the primed fields E0 ðrÞ and H 0 ðrÞ are expressed as E0 ðrÞ ¼ ½e0x ðz; j; wÞux þ e0y ðz; j; wÞuy þ e0z ðz; j; wÞuz exp½ijðx cos w þ y sin wÞ;
ð8Þ
H 0 ðrÞ ¼ ½h0x ðz; j; wÞux þ h0y ðz; j; wÞuy þ h0z ðz; j; wÞuz exp½ijðx cos w þ y sin wÞ;
ð9Þ
where w is an angle in the xy plane relative to the x axis, and j is a transverse wavenumber. Eq. (6) separate into two algebraic equations and four ordinary differential equations (ODEs), when (7)– (9) are substituted therein [12]. The two algebraic equations yield [1] j ðb a Þ 0 ½h0 ðz; j; wÞ sin n þ h0y ðz; j; wÞ cos n ex ðz; j; wÞ cos v sin v; x0 a s x a s j 0 h0z ðz; j; wÞ ¼ ½e ðz; j; wÞ sin n þ e0y ðz; j; wÞ cos n; xl0 x
e0z ðz; j; wÞ ¼
ð10Þ ð11Þ
where n¼
pz w; X
s ¼ cos2 v þ ðb =a Þ sin2 v:
ð12Þ
Thereafter, the four coupled ODEs can be represented compactly as the matrix ODE d 0 iX ½f ðf; j; wÞ ¼ ½P 0 ðf; j; wÞ ½f 0 ðf; j; wÞ; df p
ð13Þ
which has been written with the normalized axial distance f
pz X
ð14Þ
372
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
as the independent variable. The electromagnetic field components are contained in the column 4vector 2 0 3 ex ðf; j; wÞ 6 e0y ðf; j; wÞ 7 7 ð15Þ ½f 0 ðf; j; wÞ ¼ 6 4 h0x ðf; j; wÞ 5: 0 hy ðf; j; wÞ The 4 4 kernel matrix in (13) is 2 3 0 xl0 0 ip X 6 ip 0 xl0 0 7 X 7 ½P 0 ðf; j; wÞ ¼ 6 4 0 5 x0 c 0 ip X x0 b ip 0 0 X s 2 A cos n sin 2v 0 B sin n cos n 6 A sin n sin 2v 0 B sin2 n þ6 2 4 C sin n cos n C cos n 0 C sin2 n C sin n cos n A sin n sin 2v
3 B cos2 n B sin n cos n 7 7; 5 0 A cos n sin 2v
ð16Þ
where A¼
jðb a Þ ; 2a s
B¼
j2 ; x0 a s
C¼
j2 : xl0
ð17Þ
3. Series solution method Following Lakhtakia and Weiglhofer [12], the exact analytic (but not necessarily closed-form) solution of (13) can be written as the column-vector polynomial series 1 X ½f 0 ðf; j; wÞ ¼ ½pm ðj; wÞfm ; ð18Þ m¼0
which converges for all f because ½P 0 ðf; j; wÞ is a holomorphic function of f. For f ¼ 0, the polynomial series yields ½p0 ðj; wÞ ¼ ½f 0 ð0; j; wÞ; while the remaining coefficient 4-vectors on the right side of (18) are computed with the recurrence relation ½pm ðj; wÞ ¼
m1 1 iX X ½P 0 ðj; wÞ½pn ðj; wÞ; m p n¼0 mn1
m ¼ 1; 2; 3; . . .
The 4 4 matrixes in this recurrence relation are as follows: 2 3 0 xl0 0 ip X 6 ip 0 xl0 0 7 X 7 ½P 00 ðj; wÞ ¼ 6 4 0 5 x0 c 0 ip X x0 b ip 0 0 X s 2 3 B sin 2w B2 ½1 þ cos 2w A sin 2v cos w 0 2 B 6 A sin 2v sin w 0 ½1 cos 2w B2 sin 2w 7 2 7; þ6 C C 4 sin 2w 5 ½1 þ cos 2w 0 0 2 2 C C2 ½1 cos 2w sin 2w A sin 2v sin w A sin 2v cos w 2
ð19Þ
ð20Þ
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
3 A sin 2v sin w 0 B2m1 cos 2w B2m1 sin 2w imþ1 6 0 B2m1 sin 2w B2m1 cos 2w 7 6 A sin 2v cos w 7; ½P 0m ðj; wÞ ¼ m1 m1 4 5 0 0 m! C2 cos 2w C2 sin 2w m1 m1 C2 sin 2w C2 cos 2w A sin 2v cos w A sin 2v sin w
373
2
m ¼ 1; 3; 5; . . . ; ð21Þ
2
A sin 2v cos w m 6 i A sin 2v sin w ½P 0m ðj; wÞ ¼ 6 m! 4 C2m1 sin 2w C2m1 cos 2w
0 0 C2m1 cos 2w C2m1 sin 2w
B2m1 sin 2w B2m1 cos 2w 0 A sin 2v sin w
3 B2m1 cos 2w B2m1 sin 2w 7 7; 5 0 A sin 2v cos w
m ¼ 2; 4; 6; . . . ð22Þ 0
For solving boundary value problems, it is convenient to express the column 4-vector ½f ðf; j; wÞ in terms of a 4 4 matrix ½M 0 ðf; j; wÞ (called the matrizant), such that [1] ½f 0 ðf; j; wÞ ¼ ½M 0 ðf; j; wÞ½f 0 ð0; j; wÞ:
ð23Þ
It then follows that ½M 0 ðf; j; wÞ ¼
1 X
½R0m ðj; wÞfm ;
ð24Þ
m¼0
with ½R00 ðj; wÞ ¼ ½I; where ½I is the 4 4 identity matrix, and m1 1 iX X ½R0m ðj; wÞ ¼ ½P 0 ðj; wÞ½R0n ðj; wÞ; m p n¼0 mn1
ð25Þ
m ¼ 1; 2; 3; . . .
ð26Þ
4. Performance of series solution method A computer program, with double precision arithmetic used throughout, was written in FORTRAN 77 to calculate the matrizant ½M 0 ðf; j; wÞ. The program was executed on a desktop PC. In contrast to our earlier investigation on CLCs [11], in which terms were added in the series (24) for ½M 0 ðf; j; wÞ until all 16 elements reached a desired level of precision, the eight planewave remittances (i.e., four reflectances and four transmittances) of chiral STF/SLCs of finite thickness in the circularly polarized basis were calculated [1] and used to study convergence in the present study. It was found that, due to the complex nature of these dielectric HBMs in contrast to CLCs, the requirement that all elements in ½M 0 ðf; j; wÞ converge to a predetermined relative precision is unrealistic. In view of the use of these materials as optical elements, the planewave remittances serve as more practical measures of the accuracy of ½M 0 ðf; j; wÞ anyhow. Other target functions, such as the Jones matrix [2, Section 4.1.1], could have been used in lieu of the remittances; but we do not expect the general conclusions obtained here to be significantly altered. 4.1. Algorithm development for large f In the Bragg regime, the planewave reflectances saturate as the thickness of the chiral STF/SLC increases [15]. We chose the thickness as 60X – i.e., we computed the matrizant for f 60p – for the numerical
374
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
results reported here. This, however, presented a problem in calculating ½M 0 ðf; j; wÞ by the series solution method. The series (24) converges too slowly for f J p to be of any practical use for CLCs [11], but that problem is easily circumvented. Because the angle of rise is zero for CLCs, the f-period of the kernel matrix ½P 0 ðf; j; wÞ is p. Any f > 0 can be expressed as f ¼ ð‘ þ qÞp, where 0 < q < 1 and ‘ ¼ 0; 1; 2; . . . The matrizant for an arbitrary value of f can then be written as the product ½M 0 ðf; j; wÞ ¼ ½M 0 ðqp; j; wÞ½M 0 ðp; j; wÞ‘ . Thus, one never needs to directly use the series to calculate ½M 0 ðf; j; wÞ for f > p. The situation for chiral STFs/SLCs is not quite so simple. The f-period of the kernel matrix is 2p for j 6¼ 0, but only p when j ¼ 0 [19]. A compact algorithm should accommodate both conditions. Now, ½M 0 ðf; j; wÞ can be easily computed by the series (24) when 0 6 f 6 p, but not when p < f 6 2p. In other words, if we rewrite ½M 0 ðp þ Df; j; wÞ ¼ ½Q0 ðDf; j; wÞ½M 0 ðp; j; wÞ;
0 < Df 6 p;
ð27Þ
the key issue is the computation of the auxiliary matrizant ½Q0 ðDf; j; wÞ. An expression for the auxiliary matrizant can be derived as follows: After a displacement of f ¼ p, the permittivity dyadic ðfÞ is rotated 180° about the z axis. If the coordinate system is also rotated by 180° about the z axis, the elements of ðfÞ take on the same numerical values in the second half-period as in the first half-period. In the rotated coordinate system however, w becomes p þ w. Therefore, ½Q0 ðDf; j; wÞ ¼ ½I ½M 0 ðDf; j; p þ wÞ ½I ¼ ½M 0 ðDf; j; p þ wÞ; 0 < Df 6 p; ð28Þ where the two factors of ð½IÞ represent a similarity transformation. Any f > 0 can be represented as f ¼ ð2j þ ‘ þ qÞp, where 0 < q < 1; either ‘ ¼ 0 or ‘ ¼ 1; and j ¼ 0; 1; 2; . . . Therefore, 8 j > < ½M 0 ðqp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ; ‘ ¼ 0; j ð29Þ ½M 0 ðf; j; wÞ ¼ > : ½M 0 ðqp; j; p þ wÞ½M 0 ðp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ; ‘ ¼ 1: The two conditional parts of the foregoing equation can be combined to yield the equation j ‘ ½M 0 ðf; j; wÞ ¼ ½M 0 ðqp; j; ‘p þ wÞ½M 0 ðp; j; wÞ ½M 0 ðp; j; p þ wÞ½M 0 ðp; j; wÞ ;
ð30Þ
which holds for both ‘ ¼ 0 and ‘ ¼ 1. Thus, the matrizant can be computed using the series (24) for any value of the normalized axial distance f. 4.2. Model system The number of terms in the series (24) that can be summed is limited to 170 on our desktop PC, due to overflow in calculating m!. Therefore, the 170-term sum was used as a best estimate of the true matrizant by which the convergence of the series (24) could be judged. To characterize the performance of the series solution method, the summation was truncated at m ¼ N , and the value Nreq of N required to yield all remittances to within 0.1% of the values predicted by the 170-term calculation of ½M 0 ðf; j; wÞ was determined for a range of parameters describing both the incident plane wave and the constitutive properties. The starting model system for the present study is delineated by: X ¼ 200 nm, a ¼ 2, b ¼ 3, c ¼ 2:4, and v ¼ 25°. The wavelength of the incident light in free space is k0 ¼ 600 nm and the plane of incidence is described by w ¼ 0°. Except for the determination of Nreq vs. f, the thickness of the chiral STF/SLC was set
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
375
equal to 61:5X. This value was chosen so that not only was there saturation in the values of the remittances with respect to f [1], but also that our algorithm was tested on all three classes of matrixes appearing on the right side of (30). The transverse wavenumber j ¼ ð2p=k0 Þ sin h, where h is the angle of planewave incidence with respect to the z axis. Axial propagation occurs inside the chiral STF/SLC for h ¼ 0°, while oblique propagation requires h 2 ð0°; 90°Þ. 4.3. Nreq vs. f Fig. 1 shows the curves of Nreq vs. f for h ¼ 0°; 30°; 60°, 89°, and 0 < f K 200. Calculation of the set of matrizant factors in (30) has an algorithmic periodicity (not to be confused with the structural periodicity) in f of 4p. However, as f increases, the last factor in (30) is raised to ever higher powers j – thereby leading to the conjecture that, in order to maintain the same level of accuracy in the final result, higher precision in this factor might be required as f is increased. This does not appear to be the case. After a very sharp rise in Nreq for small f, the curves almost level off. A small amount of oscillation is observed in every curve. It is most pronounced at h ¼ 60° and amounts to approximately a 10% fluctuation. The maximum value of Nreq for all h chosen is slightly greater than 65. Parenthetically, the noisiness of data in Fig. 1 (as also Figs. 2–5 and 8–11) can be attributed to two reasons: (i) the periodic nonhomogeneity of HBMs, which is responsible for the differential transmissivity characteristic of the circular Bragg phenomenon; and (ii) because the quantities plotted either vary or were allowed to vary in discrete steps. 4.4. Nreq vs. X Fig. 2 shows the curves of Nreq vs. X 2 ½50; 600 nm, for h ¼ 0°; 30°; 60°, and 89°. The curves let us conclude that Nreq X, which can be explained in the following manner: As the period 2X increases, the morphological length-scale becomes comparable to the wavelength, thereby enhancing the complexity of the spatial distribution of the electromagnetic field. The increased complexity entails a larger number of
70
60
Nreq
50
40
30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10
0 0
50
100
150
200
ζ Fig. 1. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the normalized axial distance f ¼ pz=X for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and v ¼ 25°.
376
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 160 θ
140
0 deg. 30 deg. 60 deg. 89 deg.
120
Nreq
100 80 60 40 20 0 0
200
400
600
Fig. 2. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the half-period X for: k0 ¼ 600 nm, w ¼ 0°, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
terms in the series (24) for convergence. Additionally, we note that Nreq for h ¼ 89° is a little more than three times that for h ¼ 0°. 4.5. Nreq vs. c Fig. 3 displays the curves of Nreq vs. c , with c ranging from considerably below a to considerably above a . These curves do not evince an easily identifiable pattern. The curves for lower values of h are
100
θ 0 deg. 30 deg. 60 deg. 89 deg.
90 80 70
Nreq
60 50 40 30 20 10 0 1.5
2
εc
2.5
3
Fig. 3. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus c for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
377
nearly flat. Those for higher values of h show a minimum in Nreq at c ¼ 2:7, and the depth of this minimum increases with the angle of incidence. As with the other sets of calculations reported here, Nreq rises with increasing h. 4.6. Nreq vs. v The relationship between Nreq and v, shown in Fig. 4, is nearly flat. There is one notable exception at v ¼ 45° where there is a very sharp peak in the h ¼ 0° curve. This peak corresponds to the pseudoisotropic point (i.e., c ¼ b =s) [20]. At the pseudoisotropic point, the Bragg regime vanishes identically for normal incidence and is seriously weakened for oblique incidence. The general trend of larger Nreq for larger h is also evident in Fig. 4. 4.7. Nreq vs. precision The value of Nreq was calculated for relative error in the remittances ranging four orders of magnitude from 105 to 101 . The results of these calculations are presented in Fig. 5. For the chosen values of h, the maximum change in Nreq over the large range of relative precision in all remittances is about 60%. 4.8. Calculation time The calculations presented here were carried out on a desktop PC with a clock speed of 1.8 GHz. The time required to calculate the matrizant was determined by timing 100 iterations of the calculation. This number of iterations gave consistent results when the calculation time Tseries (in millisecond) was remeasured, but was not so long as to be tiresome. The time to complete the calculation depends only on the number of terms calculated in the series. The variation of Tseries with respect to Nreq is displayed in Fig. 6. The data were fit with a quadratic with excellent results; indeed, on the scale of the figure, the collected data and the quadratic fit
80 70 60
Nreq
50 40 30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10 0 0
50
χ (deg.)
100
Fig. 4. Series solution method. Number of terms required to give a relative precision of 0:1% in remittances versus the rise angle v for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and z ¼ 61:5X.
378
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 80 70 60
Nreq
50 40 30 θ
20
0 deg. 30 deg. 60 deg. 89 deg.
10 0
10-5
10-4
10-3 Precision
10-2
1
10-1
Fig. 5. Series solution method. Number of terms required for various values of relative precision in remittances for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
250
Tseries (ms)
200
150
100
50
0 0
50
100
150
200
Nreq Fig. 6. Series solution method. Calculation time versus the number of terms.
2 Tseries ¼ 8 103 Nreq þ 7:6 102 Nreq
ð31Þ
are indistinguishable. Considering the way in which the matrix coefficients of the series (24) are calculated recursively as per (26), we might have anticipated the quadratic time-dependence of the series solution method. For the coefficient matrix of the mth term in the series, m matrix multiplications are required on the right side of (26). A total of ð1=2ÞNreq ðNreq þ 1Þ matrix multiplications must be carried out for a series with 2 Nreq terms. For large Nreq , the number of matrix multiplications is thus proportional to Nreq . Except for very large values of X, the number of terms required for all calculations reported here was Nreq < 90. Therefore, the calculation time by the series solution method was Tseries 6 72 ms.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
379
5. Piecewise homogeneity approximation method The piecewise homogeneity approximation method has been successfully applied to chiral STFs/SLCs [1,14]. In this method, the material – which is nonhomogeneous along the z axis – is divided into a series of slices (or pieces) perpendicular to that axis, each of thickness Dzsl . Each slice is treated as if it were homogenous with its permittivity dyadic taken to be that at the center of the slice. In a medium with uniform permittivity dyadic, ½P 0 ðf; j; XÞ becomes independent of f [9], and the solution to (13) can be written in closed form as [21, Section 2.5] iX 0 ½P ðj; wÞDf ½f ðf; j; wÞ; Df ¼ f0 f: ½f 0 ðf0 ; j; wÞ ¼ exp ð32Þ p Therefore, if fj1 represents the normalized axial distance of the start of the jth slice, and Dfsl is the normalized slice thickness, the transfer of electromagnetic fields across this slice is expressed by iX 0 Df ½P fj1 þ sl ; j; w Dfsl ½f ðfj1 ; j; wÞ; j ¼ 1; 2; 3; . . . ½f 0 ðfj ; j; wÞ ¼ exp ð33Þ p 2 Denoting the matrizant describing propagation through the jth slice by iX 0 Df ½M 0j ðj; wÞ ¼ exp ½P fj1 þ sl ; j; w Dfsl ; p 2
ð34Þ
we see that propagation through a normalized axial distance f ¼ Nsl Dfsl is described by ½M 0 ðf; j; wÞ ½M 0N ðj; wÞ½M 0N sl
sl 1
ðj; wÞ ½M 01 ðj; wÞ;
ð35Þ
where the number of slices Nsl must be adequate for a convergent solution. In order to execute the calculation of ½M 0 ðf; j; wÞ via (35), the matrix exponential in (34) must be determined. Now, as n 1 X iX 0 iX n ½P ðfj1 þ Dfsl =2; j; wÞDfsl ¼ ½I þ Dfsl ½P 0 ðfj1 þ Dfsl =2; j; wÞ exp ð36Þ p p n¼1 is an infinite series, the sum matrix on the right side of (36) was computed term by term. Following Venugopal and Lakhtakia [1], the summation was stopped at that value of n for which the changes in all elements of the sum matrix were less than or equal to one part in 1012 . Parenthetically, we note that the left side of (36) can be evaluated via diagonalization of ½P 0 ðfj1 þ Dfsl =2; j; wÞ, following eigenvalue analysis [22], provided no eigenvalue has a geometric multiplicity other than unity [21, Section 2.5].
6. Performance of piecewise homogeneity approximation method The accuracy of the piecewise homogeneity approximation method is controlled by Dfsl . However, there is no general way of knowing, in advance, how thin to make the slices in order to obtain a particular accuracy in the calculation. To obtain a desired degree of accuracy, one must repeat the calculation several times, increasing the number of slices each time, and watch for stabilization of the remittances. We carried out a series of calculations to look at how the required Dfsl varies as a function of the parameters describing the incident plane wave and the chiral STF/SLC. The accuracy of the approximation was judged by comparing results to those of the 170-term series solution of Section 4: all remittances had to be within 0.1% of those from the series solution. The parameters and the range of independent variables studied are the same as those for studying the series solution method.
380
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
6.1. Algorithm development Although there is no difficulty in calculating the entire matrizant directly using the piecewise homogeneity approximation method for any f, the calculation time for large f can be reduced substantially by exploiting the periodicity of chiral STFs/SLCs. The factorization is simpler than for the series solution method, since the piecewise computation evinces no problems for ½M 0 ð2p; j; wÞ. Any f > 0 can be expressed as f ¼ 2ð‘ þ qÞp, where 0 < q < 1 and ‘ ¼ 0; 1; 2; . . . The matrizant for an arbitrary value of f can then be obtained as the product ‘
½M 0 ðf; j; wÞ ¼ ½M 0 ð2qp; j; wÞ½M 0 ð2p; j; wÞ :
ð37Þ
In order to determine Dfsl , Nsl was increased until the desired accuracy in the remittances was obtained. The number of slices was partitioned between the two intervals 0 6 f 6 2qp and 0 6 f 6 2p – associated with the two factors on the right side of (37) – in proportion to the size of each interval and rounded to the nearest integer. The average slice thickness was then taken to be hDzsl i ¼ ð1 þ qÞ
2X ; Nsl
ð38Þ
which quantityÕs variation with various incidence and constitutive parameters was examined. 6.2. hDzsl i vs. h Fig. 7 shows hDzsl i vs. h for the starting model system of Section 4.2. Most notable is the very large spike at h ¼ 0°. This is to be expected since it is possible to calculate ½M 0 ðf; j; wÞ for h ¼ 0° with just one slice for any f [15]. As we chose f ¼ 61:5p and our algorithm is based on (37), at least two slices are required; accordingly, hDzsl i ¼ 350 nm for X ¼ 200 nm, just as found. Because hDzsl i for h ¼ 0° is about two orders of magnitude greater than for h only slightly greater than 0° for all studies reported here, the curves for h ¼ 0° are not shown in the subsequent figures in order to not obscure the results for other values of h.
400 350
<∆Zsl> (nm)
300 250 200 150 100 50 0 0
10
θ (deg.)
20
Fig. 7. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus angle of incidence h for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
381
6.3. hDzsl i vs. f Fig. 8(a) and (b) shows the variation of hDzsl i with f for h ¼ 30°, 60° and 89°. The h ¼ 30° curve in Fig. 8(a) is essentially flat. The h ¼ 89° curve shows an initial steep decline tapering off to a much more gradual decline as a function f, while the curve for h ¼ 60° in Fig. 8(b) shows a nearly steady decline. This greater effort required to calculate the matrizant (smaller slice thickness) as f increases, is contrary to what was observed for the series solution with Nreq relatively constant. Furthermore, all three curves have a significant amount of fluctuation about the general trend, with the h ¼ 60° curve showing the largest amount. hDzsl i changes by approximately a factor of 5 in one cycle in the most extreme fluctuation. This is more than an order of magnitude greater fluctuation than the 10% effect seen in Nreq in the series solution. Over the entire range of calculations, hDzsl i changes by greater than a factor of 20. In contrast, the change in Nreq for the series solution, excluding a few points at f < 1, changes only by roughly a factor of 3 over the entire range of calculation. These observations allow us to conclude that the error in replacing a thin nonhomogeneous slice by a homogeneous slice in the piecewise homogeneity approximation method accumulates as the normalized axial thickness f increases. In order to offset the accumulation of error, hDzsl i has to be reduced. 6.4. hDzsl i vs. X The variation of hDzsl i with X is shown in Fig. 9 for h ¼ 30°, 60°, 89°. For h ¼ 60° and 89°, there is a general tendency to larger hDzsl i as X increases. The h ¼ 30° curve shows a similar trend but with a broad dip near X ¼ 500 nm, the dip corresponding to virtually null values of the cross-polarized transmittances and the co-polarized reflectances. Leaving aside that broad dip, we conclude that the number of slices increases with X. This is in line with the series solution method for which Nreq (and hence the computational effort) increases with X for all angles of incidence. For both methods, the computational effort is more at higher h.
7
7 θ
6
6 θ
30 deg.
4 3
4 3
2
2
1
1 0
0 0
(a)
60 deg.
5
89 deg.
<∆Zsl> (nm)
<∆Zsl> (nm)
5
50
100
ζ
150
0
200
(b)
50
100
ζ
150
200
Fig. 8. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the normalized axial distance f ¼ pz=X for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°.
382
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386 10 9
θ
8
30 deg. 60 deg. 89 deg.
<∆Zsl> (nm)
7 6 5 4 3 2 1 0 0
200
400
600
Ω (nm)
Fig. 9. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the half-period X for: k0 ¼ 600 nm, w ¼ 0°, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
6.5. hDzsl i vs. c There is no general trend evinced by the hDzsl i vs. c curves in Fig. 10. Oscillations in these curves, however, are quite dramatic – most prominently for h ¼ 30° when c K 2:25. The oscillations are in marked contrast to the dependence of Nreq on c in Fig. 3 for the series solution method. The oscillations in the hDzsl i vs. c curve coincide with oscillations in the remittances inside the Bragg regime. Furthermore, the drop for h ¼ 30° when c J 2:25 corresponds to drops in i(i) the co-polarized reflectance when the incident plane wave is right circularly polarized, and (ii) the co-polarized transmittance when the incident plane wave is left circularly polarized, 9 θ
8
30 deg. 60 deg. 89 deg.
7
<∆Zsl> (nm)
6 5 4 3 2 1 0 1.5
2.0
εc
2.5
3.0
Fig. 10. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus c for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, v ¼ 25°, and z ¼ 61:5X.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
383
outside the Bragg regime. The correspondence between the oscillations in Fig. 10 and the remittances is less striking when h ¼ 60°, and virtually vanishes when h ¼ 89°. The hDzsl i-range for calculations at all h is quite large. The smallest value of hDzsl i in the h ¼ 89° curve is 0.0772 nm, while the largest value in the h ¼ 30° curve is 8.14 nm. The two values differ by slightly more than two order of magnitude. This is very much greater than the range of values seen in Nreq of the series solution for which the maximum value of Nreq differs from the minimum only by a factor of roughly 4. 6.6. hDzsl i vs. v The curves of hDzsl i vs. v are shown in Fig. 11 for h ¼ 30°, 60°, and 89°. Again, the fluctuations in these curves are extremely pronounced. Some, but not all, features of the hDzsl i vs. v curves can be correlated to remittance variations. Recall that, for the series method, Nreq vs. v has only one peak at v ¼ 45° (the pseudoisotropic point) when h ¼ 0°. Over the range of calculations illustrated in Figs. 4 and 11, hDzsl i changes by a factor of nearly 47, while the maximum value of Nreq differs from the minimum by less than a factor of 4. 6.7. hDzsl i vs. precision Fig. 12(a) shows the relationship between hDzsl i and relative precision in the remittances which range from 101 to 105 . Unlike the series solution method – for which Nreq varies very little with relative precision (see Fig. 5) – hDzsl i depends very strongly on relative precision. The log–log plot in Fig. 12(b) shows that hDzsl i varies as the square of the relative precision. From a relative precision of 105 to 101 , hDzsl i changes by approximately two orders of magnitude; in contrast, the maximum change in Nreq is by 60% in Fig. 5. 6.8. Calculation time Calculation times with the piecewise homogeneity approximation method were determined in the same way as for the series solution method. There are two main factors which influence the calculation time Tpiece 8 θ
7
30 deg. 60 deg. 89 deg.
<∆zsl> (nm)
6 5 4 3 2 1 0 0
20
40
χ (deg.)
608
0
Fig. 11. Piecewise homogeneity approximation method. Slice thickness required to give a relative precision of 0:1% in remittances versus the rise angle v for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, and z ¼ 61:5X.
384
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
30 10 1
20
<∆Zsl> (nm)
<∆Zsl> (nm)
θ 30 deg. 60 deg. 89 deg.
10 0 θ
10
0 -5 10
(a)
30 deg. 60 deg. 89 deg.
10-1
10- 4
10- 3
10- 2
10-2
10- 1
10-5
10-4
(b)
Precision
10-3
10-2
10-1
Precision
Fig. 12. Piecewise homogeneity approximation method. Slice thickness required for various values of relative precision in remittances for: k0 ¼ 600 nm, w ¼ 0°, X ¼ 200 nm, a ¼ 2:0, b ¼ 3:0, c ¼ 2:4, v ¼ 25°, and z ¼ 61:5X.
for the piecewise method: X and hDzsl i. Measurements of Tpiece were made for hDzsl i ¼ 0:1, 0.5, 1.0, 5.0, and 10.0 nm, with X ranging from 50 nm to 600 nm. The results for hDzsl i ¼ 0:1, 1.0, and 5.0 nm are reported in Fig. 13. As can be deduced from this figure, Tpiece varies linearly with X for all values of hDzsl i. The dependences of the slopes of the Tpiece vs. X curves on hDzsl i can be represented by a simple power law. A reasonable approximation for calculation time emerges from Fig. 13 as Tpiece ¼ 0:20XhDzsl i0:89 ;
ð39Þ
where Tpiece is in millisecond and hDzsl i in nm. This formula predicts the calculation time within 10% for most cases, although the error may be as much as 25% at low values of X and high values of hDzsl i.
1000 <∆ Zsl> 0.1 nm 1.0 nm 5.0 nm
800
Tpiece(ms)
600
400
200
0 0
200
400
Ω (nm)
600
Fig. 13. Piecewise homogeneity approximation method. Calculation time versus X for various values of hDzsl i.
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
385
7. Performance comparison The series solution method and the piecewise uniform permittivity method offer two viable means of determining the propagation of light in chiral STFs and chiral SLCs. The former method is exact and analytic, whereas the latter method is heuristic. For k0 > 3X and the chosen values of the constitutive parameters, planewave remittances can be calculated to a precision of 0:1% with less than 90 terms in the series (24). With this number of terms, calculations can easily be performed on a desktop computer, and a single calculation of ½M 0 ðf; j; wÞ appears instantaneous to the user with a calculation time 672 ms. The number of terms required for higher precision is only modestly higher, still well within the range accessible by desktop computers. To increase the relative precision from 0:1% to 0:001%, the number of terms in (24) only needs to be increased by approximately 30% in the worst case that we studied, resulting in a 70% increase in calculation time. The required number of terms Nreq for a given accuracy is very weakly dependant on the parameters describing the material and the incident plane wave. A general trend is that Nreq increases with angle of incidence h. These conclusions are true both inside as well as outside the Bragg regime. With the exception of the h ¼ 89°, all Nreq vs. f curves in Fig. 1 pass through the Bragg regime. There is no substantial change evident in Nreq as this occurs. In Figs. 2–5, the h ¼ 0° and h ¼ 30° curves are in the Bragg regime, while the curves for other values of h are not. No qualitative difference is apparent between those curves that are in the Bragg regime, and those that are not. The performance of the piecewise homogeneous permittivity approximation method appears to be much more variable. The slice thickness required to obtain 0.1% precision in remittances fluctuates dramatically. For all calculations made with 0.1% precision and X ¼ 200 nm, hDzsl i varies between 0.0772 nm and 8.14 nm. This implies that the calculation time varies roughly between 6 ms and 400 ms. To increase relative precision form 0.1% to 0.001%, hDzsl i must be reduced about an order of magnitude, increasing the calculation time by nearly a factor of 8. The calculation times with both the series solution method and the piecewise homogeneity approximation method vary with the parameters describing the helicoidal bianisotropic medium and the incident plane wave. A general conclusion is that calculation times for the former method vary less erratically than for the latter method. For the following discussion, all calculation times were obtained using (31) and (39). The calculation time with either method increases with the angle of incidence of the plane wave. With the piecewise homogeneity approximation method, there is a large change in the first few degrees away from normal incidence, amounting to nearly an order of magnitude. For our starting model system, the piecewise homogeneity approximation method is roughly 10 times faster than the series solution method at h ¼ 0°; the two methods are comparable for h J 2°; but the piecewise homogeneity approximation method performs at about half the speed of the series solution method at h ¼ 89°. The calculation time with either method increases with the structural period. The rate of increase of the calculation time with the structural period is larger for the series solution method than for the piecewise approximation method. For the model system, the calculation times for the two methods are very similar when X ¼ 50 nm, but the piecewise homogeneity approximation method is faster (except at h ¼ 89°) when X ¼ 600 nm. Thus, for h ¼ 60°, Tpiece ¼ 12 ms at X ¼ 50 nm and 45 ms at X ¼ 600 nm; whereas Tseries ¼ 9 and 120 ms, respectively. Over the range of X studied, neither method appears to have a clear advantage over the other. The calculation times for the two methods clearly differ as functions of f. Except at extremely small f, Tseries is independent of f. This is true of the piecewise homogeneity approximation method at h ¼ 30° as well, but at larger angles Tpiece increases with f. At h ¼ 89°, for instance, from f ¼ 5 to f ¼ 190, Tpiece changes from about 20 to 100 ms. The series solution method, on the other hand, has a roughly constant calculation time of about 35 ms.
386
J.A. Polo Jr, A. Lakhtakia / Optics Communications 230 (2004) 369–386
No general trends regarding the calculation times emerged with respect to either the ratio c =a or v, except that the larger fluctuations in hDzsl i than in Nreq imply that Tpiece varies more erratically than Tseries with respect to these parameters.
8. Concluding remarks The exact and analytic series solution method can be applied only if the 4 4 kernel matrix in (13) is a holomorphic function of f, while the heuristic piecewise homogeneity approximation method is definitely more versatile. However, as we have shown here, the former method does offer some advantages when determining the oblique propagation characteristics of light in chiral STFs and chiral SLCs. The number of terms required in the series (24) to calculate planewave remittances to a given precision remains relatively constant over a wide range of parameters. The same cannot be said of the average slice thickness needed to implement the piecewise homogeneity approximation method. Also, to calculate the remittances for a specific set of geometric and constitutive parameters to a predetermined level of accuracy, the entire piecewise calculation has to be repeated several times with different values of Nsl . Thus, the actual calculation time for a single calculation with the piecewise homogeneity approximation method could be significantly larger than those reported here for which the slice thicknesses were pre-determined. In contrast, the convergence of the series solution method can be monitored as the series is summed and the process terminated when the desired accuracy has been achieved. This could be a significant advantage when a large set of calculations is to be performed.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
V.C. Venugopal, A. Lakhtakia, Proc. R. Soc. Lond. A 456 (2000) 125. S. Chandrasekhar, Liquid Crystals, second ed., Cambridge University Press, Cambridge, UK, 1992. P.G. de Gennes, J. Prost, The Physics of Liquid Crystals, second ed., Clarendon Press, Oxford, UK, 1993. S.D. Jacobs (Ed.), Selected Papers on Liquid Crystals for Optics, SPIE Press, Bellingham, WA, 1992. I.J. Hodgkinson, Q.h. Wu, Adv. Mater. 13 (2001) 889. A. Lakhtakia, Mater. Sci. Eng. C 19 (2002) 427. E.I. Kats, Sov. Phys. JETP 32 (1971) 1004. R. Nityananda, Mol. Cryst. Liq. Cryst. 21 (1973) 315. A. Lakhtakia, W.S. Weiglhofer, Microw. Opt. Technol. Lett. 6 (1993) 804. A. Lakhtakia, W.S. Weiglhofer, Microw. Opt. Technol. Lett. 12 (1996) 245. J.A. Polo Jr., A. Lakhtakia, Microw. Opt. Technol. Lett. 35 (2002) 397. A. Lakhtakia, W.S. Weiglhofer, Proc. R. Soc. Lond. A 453 (1997) 93, corrections: 454 (1998) 3275. M. Schubert, C.M. Herzinger, Phys. Stat. Sol. (a) 188 (2001) 1563. I. Abdulhalim, R. Weil, L. Benguigui, Liq. Cryst. 1 (1986) 155. V.C. Venugopal, A. Lakhtakia, in: O.N. Singh, A. Lakhtakia (Eds.), Electromagnetic Fields in Unconventional Materials and Structures, Wiley, New York, NY, 2000 (Chapter 5). J.A. Sherwin, A. Lakhtakia, Opt. Commun. 222 (2003) 305. A. Lakhtakia, W.S. Weiglhofer, IEE Proc.-Microw. Antennas Propag. 144 (1997) 57. C.W. Oseen, J. Chem. Soc. Faraday Trans. II 29 (1933) 883. € bertrag. 53 (1999) 287. A. Lakhtakia, V.C. Venugopal, Arch. Elektron. U A. Lakhtakia, Microw. Opt. Technol. Lett. 34 (2002) 367. H. Hochstadt, Differential Equations: A Modern Approach, Dover Press, New York, NY, 1975 (Section 2.5). M. Schubert, Phys. Rev. B 53 (1996) 4265.