Waveguide Characteristics Near the Second Bragg Condition

We show that in an optical waveguide with no material losses at wavelengths near the second-order Bragg condition, there exists two pairs of modes. One pair has identical propagation constants but have different attenuation coefficients, while a second pair with identical propagation constants (different from the first pair) has different attenuation coefficients. The four attenuation coefficients may have either a positive value, representing power leaking out of a waveguide mode or a negative value, representing power from an external source leaking into a mode. Moreover, a mode with a positive (negative) attenuation before the second Bragg condition, has a negative (positive) attenuation after the second Bragg condition. Radiation near the second-order Bragg condition of a periodic waveguide typically occurs at an angle perpendicular or nearly perpendicular to the propagation direction of the waveguide because the scattering centers have a period equal to or close to the period of the longitudinal propagation constant of the mode. In this paper, stable numerical solutions for the modes of periodic dielectric structures are developed using Floquet-Bloch theory. One primary focus in this paper illustrates a unique method of analyzing such modes using eigenvectors forming a Hilbert space, allowing for expansion of arbitrary vectors and their derivatives used for calculations such as that involving group velocities. The dimension of the vector space is determined by the number of space harmonics used in the solution of the Floquet-Bloch equations. The accuracy of the numerical solutions is affected by the number of space harmonics; as that number is increased, the number of waveguide partitions must be increased to maintain a given accuracy.

teristics that are considerably different from similarly con- 23 structed conventional waveguides. A mode propagating, say 24 in the positive z direction, exp (jωt − γ z), with periodicity 25 along z, is identified by its complex propagation constant, 26 γ = α + jβ, where β = 2π/λ g , α is the attenuation, 27 λ g is the mode wavelength in the waveguide, and ω is the 28 real angular frequency of the mode. The attenuation α is 29 The associate editor coordinating the review of this manuscript and approving it for publication was Derek Abbott . related to the loss (or gain) of mode power due to: 1) power 30 coupled out of (or into) the periodic waveguide by leaky 31 modes; 2) power transfer between co-directional and contra-32 directional propagating modes, depending on grating period 33 [1]; and 3) material loss or gain. In this paper we assume that 34 all material losses and gains are zero. 35 The Floquet-Bloch method provides a prescription for ana- 36 lyzing modes in periodic waveguides that extend to infinity 37 in the direction of propagation. A key point in this paper, 38 as discussed in Section III-F, is that computation of the 39 modal attenuation coefficient, α, can only be determined by 40 specifying boundary conditions at infinite lateral distances 41 VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ that are perpendicular to the propagation direction. Boundary second Bragg condition. In particular, the matrix resulting 95 from boundary-condition settings is singular at a specific 96 real frequency or wavelength and our algorithm using the 97 Simpson method to calculate the real frequency is rigorous 98 and numerically stable for hundreds of space harmonics.

99
A somewhat different focus is the understanding of scatter-100 ing of light shinning on planar structures with etched gratings 101 on the surface or gratings embedded within the structure [4], 102 [5]. The electromagnetic field is expanded in a Fourier series 103 and a matrix is designed to account for all scattered light due 104 to an incident plane-wave with a specific field polarization, 105 and it results from the matching of Fourier components at 106 boundaries. If the matrix is non-singular, the scattered fields 107 may be determined. If the matrix is singular at the specified 108 frequency, a resonant condition, the scattered fields cannot be 109 determined at that frequency. However, because the scattering 110 matrix is obtained by matching fields at all boundaries, the 111 singularity of the matrix implies the existence of a waveguide 112 mode [6], [7]. 113 A scattering-matrix approach has also been used to address 114 appearances of leaky waveguide modes in dielectric materials 115 below a two-dimensional periodic pattern [8] and to analyze 116 crystalline slabs with asymmetric and symmetric dielectric 117 unit cells [9]. 118 Although this work has a narrow focus on the wavelength 119 range near the second Bragg condition where the mode leaks 120 normal or near normal to the propagation direction, leaky 121 modes at high-order Bragg conditions leak at different angles 122 [10], [11]. The leaky mode radiation direction is determined 123 by interaction of the space harmonics in the Floquet-Bloch 124 expansion. Assuming the periodic structure can be modeled 125 with a set of scattering centers such as shown Fig. 2, the field 126 phases of the scattering centers fix the radiation direction. 127 When phases are identical, such as the phases at the second 128 Bragg condition, the main radiation lobe is in the normal 129 direction. At higher-order Bragg conditions, the equivalent 130 scattering-center spacing increases producing a radiation pat-131 tern with grating lobes. There is a broadside beam, but grating 132 lobes appear at an angle θ, relative to the normal of the array 133 axis, which satisfies sin θ = ±λ/ . For example, a forth-134 order Bragg condition, = 2λ, has two grating lobes at 135 θ = ±π/6. Nevertheless, unlike an antenna array that may 136 radiate at broadside, a periodic dielectric waveguide does 137 not radiate at the exact second Bragg condition because the 138 amplitude of the -1 space harmonic producing the leaky wave 139 vanishes at the exact second Bragg condition. 140 Complex coupling coefficients occur if the waveguide peri-141 odicity is a result of both index and loss (or gain) periodic-142 ity [2]. Coupled mode theory (CMT) can address coupling 143 between a guided wave and radiating partial waves [12]. Dis-144 tributed feedback (DFB) grating-coupled surface-emitting 145 (GSE) lasers were characterized from an expanded CMT to 146 active second-order gratings [13]. Modeling of second-order 147 structures was used to design GSE lasers [14]. A Greens 148 function method was proposed to obtain the complex cou-149 pling coefficient of gain or loss in second-and high-order 150 gratings [15]. The Greens function method was also used to 151 analyze and design DFB lasers with second-order gratings 152 [16]. A volume current method (VCM) was used to calculate 153 FIGURE 1. Cross section of multi-layer waveguide with a single grating layer with a periodic dielectric constant. In many cases the periodic structure is formed by etching the material in Layer 4 and filling the etched trenches with the material of Layer 2. The tooth height is t h , the tooth width is t w , and the grating period is .
the radiation scattered from gratings by using a Greens func- conductor laser [23]. A thin liner layer along with a high 185 index cover layer over a simple grating can theoretically 186 reduce the length of an out-coupler grating by an order of 187 magnitude [24]. The Fourier eigenmode expansion method 188 with Floquet theory (Ref. [25]) was applied to analyze the 189 structures in Refs.

192
The method developed in this paper analyzes character-193 istics of waveguides commonly used in the integration of 194 semiconductor lasers and optical devices [24]. Dielectric 195 waveguides containing one or more special layers whose 196 dielectric constant varies with z (the propagation direction) 197 are typically fabricated with additional parallel multiple uni-198 form layers as shown in Fig. 1. As a result, the shape of 199 the mode propagating in the z direction is dependent upon 200 both x and z. In this paper, the characteristics of waveguide 201 modes in a multi-layer structure, with special layers whose 202 dielectric constant varies along the propagation direction, 203 will be developed. All layers, except the grating layers, are 204 uniform and isotropic.

206
A coordinate system is chosen such that the grating layers 207 are referenced at the transverse location x = 0 such as 208 illustrated in Fig. 1 for a simple grating formed in Layer 3. 209 The propagation direction z has an origin at the center of 210 the tooth so that the dielectric constant is symmetric about 211 the tooth center. The periodic grating of Layer 3 has two 212 regions with two different dielectric materials. The dielec-213 tric constant of the lth layer is given by κ l . Although the 214 dielectric materials of the grating layer may be formed with 215 arbitrary materials, it is typically formed with the dielectric 216 materials of the neighboring layers. For example, the central 217 tooth region usually would have a dielectric constant κ t = 218 κ 4 and the fill region would have a relative dielectric constant 219 κ f = κ 2 .

220
In the grating layer the dielectric constant is independent 221 of y and x but it is periodic with respect to z, the propagation 222 direction, and may be expanded in a spatial Fourier series 223 (The grating layer in Fig. 1 (2) 228 A non-grating layer has only the m = 0 Fourier coefficient. 229 When the grating layer thickness, t h is small compared to 230 the width of the transverse field, the periodic layer will have 231 only a slight influence (perturbation) on the shapes of the 232 transverse fields. A first-order analysis of the electromagnetic 233 fields can be obtained by replacing the periodic layer with a 234 material whose dielectric constant is the average value In the absence of material dispersion, 239 all Fourier coefficients are independent of the grating period, 240 , and the free-space wavelength, λ.

241
The Fourier expansion of the relative dielectric constant in 242 (1) can be used to extract information regarding the interac-243 tion with an electromagnetic mode that impinges on a region 244 VOLUME 10, 2022 FIGURE 2. Scattering from each of the scattering centers. When the round trip optical path length between 2 and 3 is a multiple of 2π back scattering is enhanced.
with an embedded grating such as that in Fig. 1 in Fig. 1, the Fourier coefficients for layer l = 3 are given as  of the coefficient K l m is calculated by multiplying the factors 288 in the table with the relative dielectric constant difference 289 κ t − κ f . The calculations for duty cycles of 2/3 and 3/4 can be 290 obtained using the coefficients in Table 1 by interchanging κ t 291 and κ f .

293
Early studies of periodic waveguides were focused on the 294 problems in antenna design and analysis, such as wave guid-295 ing, scattering and radiation patterns [27], [28], [29], [30]. 296 In a dielectric slab waveguide without loss or gain, in any 297 layer, there is no power loss for a bounded propagating mode. 298 Waveguide power loss can occur if there are waveguide dis-299 continuities. Since a grating tooth is a discontinuity, a peri-300 odic sequence of teeth produces significant scattering.

301
Wave phenomena such as that of electromagnetic propaga-302 tion in periodic waveguides can be described by the superpo-303 sition of all forward and backward waves which are excited 304 by the periodic structure. According to Floquet's theorem, 305 the fields can be characterized mathematically in terms of the 306 phase difference between the propagating fields and the fields 307 are discrete and are equal to the multiple of the unit reciprocal 308 grating period, K = 2π/ .

309
In general, the field distribution of a wave propagating 310 in the positive z direction in a waveguide that is periodic 311 along the propagation direction z can be expressed in terms 312 of Floquet-Bloch analysis as where R(x, z) = R(x, z + ) is a periodic function along z 315 and γ = α + jβ is the complex propagation constant. The 316 attenuation coefficient α may have a positive or a negative 317 value, however, the forward wave requires β > 0.

318
The Fourier series for R(x, z) is written as 320 which produces the field solution where γ n = γ + jnK .

323
In general, the space harmonics ψ n (x) and ψ * m are not 324 orthogonal, i.e.,

358
Substituting the first two equations into the third gives the 359 wave equation representing the field component E y .

360
The Helmholtz equation for the isotropic grating layer can 361 be written as

363
where k is the free-space wavenumber. The resulting field 364 solutions are 366 367

368
The periodic relative dielectric constant is separated into

390
A set of differential equations for the space harmonics 391 in the lth layer can be obtained by substituting (8) into the 392 Helmholtz equation, (13) and using the expansion, (1) The summation in (19) represents the discrete convolution 395 of the space harmonics with Fourier coefficients of the rela-396 tive dielectric constant.

397
If the largest coefficient, K l 0 , of the lth layer is removed 398 from the summation in (19), one gets The summation which omits the m = n term, acts as a 402 perturbation on the set of uncoupled differential equations; it 403 represents the convolution of the space harmonics with the 404 Fourier components of κ .

405
In Fig. 1, the grating tooth is rectangular, rendering κ l 406 and K l n−m independent of x. When the grating tooth is not 407 rectangular, then κ l and the corresponding K l n−m Fourier 408 coefficient have transverse x dependence, which results in a 409 set of differential equations with non-constant (with respect 410 to x) coefficients. (While differential equations with con-411 stant coefficients have well-known solutions, solutions of 412 differential equations with variable coefficients may not have 413 simple closed-form answers and one must rely on the use of 414 numerical methods.) 415 The fill regions of the grating may consist of multiple 416 material compositions. For example, Fig. 3 illustrates a single 417 grating period that consists of a tooth region and several fill 418 VOLUME 10, 2022  regions. Furthermore, the dielectric constant in the grating 419 layer is also independent of the transverse x dimension.

420
As the tooth and fill regions increase in number, the vari-

455
A main contribution of this work is identifying and using the 456 eigenvectors of the eigenvalue problem stated in (25). The 457 vectors allow computation of derivatives required for finding 458 the propagation constant via Newton's method, group veloci-459 ties and other parameters requiring the necessary derivatives. 460 One technique of solving for the fields in the grating is 461 to assume that only a finite number of space harmonics 462 in (19) play a significant role in the determination of the 463 mode. Assuming that the smallest indexed space harmonic 464 is ψ L while the largest is ψ H (the total number of space 465 harmonics is M = H − L + 1), then the infinite set of 466 differential equations reduces to a finite set of M second-467 order differential equations that may be written in vector form 468 as (For the sake of convenience, the layer subscript has be 469 dropped.) where is a diagonal matrix of γ n 's, and K is a Toeplitz 472 matrix whose entries are obtained from the set of Fourier 473 coefficients {K n } in (1). The resulting matrix P has elements 474 where δ mn is the Kronecker delta. The vector ψ is and has entries made from the individual space harmonics. 478 It should be noted that the set of M second-order equations 479 may be reduced to a set of 2M first-order equations. The set of 480 2M equations will produce a set of 2M x-dependent solutions. 481 When the matrix P is independent of x, the solution of 482 the vector differential equation, (21), can be written in terms 483 of a vector of constants, v, and a transverse propagation 484 constants, h, as The solution of (19) reduces to the algebraic eigenvalue 489 problem where σ n = h 2 n is the nth eigenvalue of P and v n 490 is the corresponding eigenvector. There are two exponential 491 functions in the solution vector, composed of exp(jhx) and 492 exp(−jhx) or a combination of the two exponential func-493 tions. The number of different pairs of lateral propagation 494 constants, (h n = ± √ σ n ), corresponds to the dimension of the 495 P matrix and is equal to M , the number of space harmonics 496 that will be used to form the solution. The set of eigenvectors 497  The mth vector v m with the corresponding eigenvalue σ m , 501 is called a right-hand eigenvector. There is a left-hand eigen-502 vector u m that has the same eigenvalue and satisfies where u H m is the transpose of the complex conjugate of u m .

505
The vectors in the set S u form a Hilbert space and are orthogo-506 nal to vectors in S v , [34] i.e., if the eigenvalues σ m = σ n , then  of v and σ with respect to γ , required when using Newton's 530 method as discussed in references [37] and [38].
The first-order results allow for the calculation of σ (1) n , 552 and v (1) n , given by where the matrix elementsp Since the set of vector equations are derived from the second-562 order differential equation, the general solution will be writ-563 ten as the sum of 2M = 2(H − L + 1) independent solutions 564 as The coefficient of v ν is the expansion variable that depends 567 only on the index ν so that a general vector, ψ in the space S v 568 will be a linear combination of the base vectors as in (27), The complete set of base eigenvectors and their derivatives 571 can be written as matrices, For multiple grating layers, say the lth layer, the matrices 574 are denoted as V l and V l γ . As opposed to non-grating layers 575 that have a single transverse wave number associated with a 576 space harmonic, in a grating layer a single transverse wave 577 number h ν does not define a space harmonic, as a space 578 harmonic is obtained from a specific combination of the set 579 of eigenvalues, {h 2 ν }. Indeed, rows of V are associated with 580 the individual space harmonics, viz.,

628
For space harmonic n, the corresponding eigenvalue, h 2 n , 629 of the lth grating layer satisfies The free-space wavenumber k and the effective index of 632 refraction of the propagating mode n e have replaced the 633 propagation constant β.

634
To investigate the solution obtained for various eigenval-635 ues, the solution of h 2 for different space harmonics can easily 636 be observed. For example, when n = 0, Because k 2 n 2 e > k 2κ l , h l 2 0 lies in the second or third 639 quadrant, depending on the value of α, which is generally 640 small compared to other terms in the expression. Splitting h l n 641 into its real and imaginary parts, h l n = h l n + jh l n , produces 642 |h l 0 | |h l 0 |, so that the corresponding field associated with 643 h will be predominately evanescent. On the other hand, near 644 the second Bragg condition, when the grating wave number 645 K ≈ k n e , the eigenvalue for the n = −1 space harmonic 646 becomes so that h l −1 ≈ k √κ l indicating the fields are oscillatory rather 649 than evanescent, corresponding to light radiating away from 650 (or into) the waveguide.

651
Employing a large number of space harmonics will of 652 course, increase the accuracy yielded by the model. On the 653 other hand, using a large number of space harmonics makes 654 the model vulnerable to numerical inaccuracies associated 655 with machine precision. (The numerical accuracy can be 656 addressed by splitting thick layers to a number of thin layers.) 657 As n becomes large, which clearly indicates the fields are evanescent. The expres-660 sion in (36) represents the field when the space harmonics are 661 coupled (the off-diagonal terms of the Toeplitz matrix K are 662 not assumed to be zero). However, the approximate solution 663 may be written as where h l n in (41) is replaced by the dominant imaginary 666 part jh l n . When the layer thickness d l makes the numerical 667 value of exp(−2h l n d l ) below machine precision, calcula-668 tions become ill-conditioned, causing numerical instability. 669 Numerical stability requires splitting thick layers to several 670 smaller layers.

671
At wave lengths near the second Bragg condition, the 672 discussion above points to the fact that all of the space har-673 monics, save n = −1, have electromagnetic fields in the 674 grating region that are predominately evanescent.

676
Exterior to the grating region, space harmonics are uncoupled 677 because the dielectric constant is independent of z, and in the 678 lth layer (19) becomes Layer-interface points x l are illustrated in Fig. 1. In the two 686 semi-infinite regions the solutions will be written specifically  wave at the 1st scattering center in Fig. 2 is 0, the phase at 730 the second scattering center will be retarded by an amount, 731 say, exp(−jϕ); the phase at the 3rd scattering center will 732 be retarded by exp(−j2ϕ). Assuming the scattering centers 733 represent a phased array, the far-field will peak at an angle of 734 −ϕ relative to the x axis, the normal to the array axis. In Fig. 5, 735 the red leaky waves occur below the second Bragg condition 736 whereas the blue leaky waves occur above the second Bragg 737 condition.   Table 2; the duty cycle is 50% and the grating period 753 VOLUME 10, 2022  (51) 789 Fig. 8(a) illustrates the real and imaginary parts of the 790 normalized group velocity of Mode I for < B and for 791 Mode II for > B . Both Mode I and Mode II leak out 792 for the grating periods shown. Fig. 8(b) shows the real and 793 imaginary parts of the normalized group velocity of Modes III 794 and IV for the specified grating periods. Both modes leak out 795 in the superstrate and leak in in the substrate for the illustrated 796 grating period range.

797
Above and below the second Bragg condition the real parts 798 of the group velocity for Modes I and II are identical, how-799 ever, their imaginary parts change sign. Namely, if Mode I and 800 Mode II regions in Fig. 8(a) are reversed, the imaginary part 801 of the group velocity will have a sign change. For example, 802 the imaginary part of the group velocity for Mode I is negative 803 for < B , whereas, the imaginary part of the group 804 velocity for Mode II is positive for the same region, < B . 805 Modes III and IV exhibit behavior similar to the behavior of 806 Modes I and II.

807
The real part of group velocity may not indicate energy 808 velocity in a periodic structure [39]. However, when a Bloch 809 mode propagates with little or no loss then {v g } approx-810 imates the energy velocity. For example, the modes of the 811 structure listed in Table 2 have attenuation coefficients that 812 are relatively small for the grating periods below the second 813 Bragg condition. The imaginary part of the group velocity 814 for the modes is approximately 0. Thus, the real part of the 815 group velocity will approximate energy velocity. On the other 816 hand, Modes III and IV have large attenuations above the 817 second Bragg condition and the corresponding real part of 818 their group velocity is large. Note that there appears to be a 819 Kramers-Kronig relationship between the real and imaginary 820 parts of the group velocity above the second Bragg condition. 821 Namely, v g peaks when v g = 0.

822
The modes are characterized as follows: the binary num-823 ber 0 represents a leak-out mode whereas the binary number 1 824 represents a leak-in mode. Thus, 00 represents a mode that 825 leaks to the superstrate as well as to the substrate. The value 826 01 represents a mode that leaks to the superstrate but leaks 827 from the substrate. Thus, Mode I is characterized by 00 before 828 the second Bragg condition, and is characterized by 11 after 829 the second Bragg condition.  clad regions will be toward the central waveguide region. 859 In particular, at z = 0, (18) impliesR(x, 0) = R * (x, 0) (for the 860 symmetric grating) so that the modes leak along x in opposite 861 directions.

862
The coefficient A l n represents the field at the base of each 863 layer while B l n represents flux or field derivative with respect 864 VOLUME 10, 2022 to x. The tangential field components at the layer interfaces 865 are E y and H z , and (14) and (16)  The results above will be applicable for all layers including Assuming the structure of Fig. 1, matching fields will 889 produce a set of nine equations for each space harmonic . . . The diagonal block matrices have band shapes (one lower 912 diagonal, and two upper diagonals) and are identical for each 913 space harmonic. For example, the B nn block for the nth space 914 harmonic in (57), as shown at the bottom of the next page.

915
For sake of clarity, a superscript is placed on the matrix ele-916 ments of V to denote that the grating is in Layer 3. If Layer 2 917 was also a grating layer, the matrix elements V 2 nn would 918 appear in columns 2 and 3. The block matrix B nn and its 919 derivative are formed from a series of sub-block matrices: 920 For the lth layer (non-grating) and nth space harmonic, the 921 modified transfer matrix T l nn and its derivative T l nn γ are For the lth grating layer and mth space harmonic, the 924 modified transfer matrices and their derivatives are The off-diagonal block matrices are composed of modified 929 transfer matrices that have six non-zero elements, which is six 930 times the number of grating layers. The off-diagonal blocks, 931 B mn , are given by The matrix 0 s is a 9×3 block of 0's while 0 r is a 9×4 block 934 of 0's. The non-zero elements line up with the V elements 935 of (57). The transfer matrix in B mn is T 3 mn . If Layer 2 is also a 936 grating layer, then 0 s would be replaced with one column of 937 0's and two columns of T 2 mn , placed at row one. The field coefficients A l n and B l n can also be partitioned or 939 blocked as follows   .4)) using a sparse matrix solver [40].

974
A node is defined as a transverse position x, where the field is 975 calculated from the solution of (60). In the derivation of (60), 976 the nodes are specified as points that lie at the interface of two 977 distinct layers, including the different layers of the grating 978 region. In general, nodes may be placed at any point in the 979 transverse direction x. 980 Fig. 10 (a) shows five node positions where the nodes are 981 placed at layer interfaces so that the node position is ξ i = x i . 982 However, Fig. 10 (b) has eight nodes with the addition of four 983 new nodes, while two nodes are moved from that shown in (a). 984 The two end nodes are anchor nodes and cannot be deleted 985 but they can be moved to cladding regions one and six. For 986 the node positions shown in Fig. 7 (b), the field quantities A 987 and B would be computed at all eight node positions for each 988 space harmonic. When a node does not lie on an interface 989 the transfer matrix will be computed as a product of two 990 single layer transfer matrices. For example, when the material 991 between x 5 and x 6 is different from the material between 992 x 4 and x 3 the transfer matrix from ξ 7 to ξ 6 would be the 993 product of two simple transfer matrices.

995
The number and location of nodes affect the accuracy of 996 the calculation of the eigenvalue γ and the corresponding 997 eigenvector (58). The characteristics of the simple dielectric 998 structure [3], [31] with an etched grating will be derived. The 999 dielectric structure described by Table 2 is a basic four-layer 1000 waveguide where the superstrate, Layer 1, is air. (This par-1001 ticular structure has been extensively investigated by many 1002 different techniques.) A grating with a 50 percent duty cycle 1003 is formed in Layer 2 where the dielectric constant, varying 1004 between 3.0 and 1.0 along the propagation direction, pro-1005 duces an average value of 2.0.
(57) VOLUME 10, 2022 FIGURE 11. Dielectric profile corresponding to Table 2  The angles for the leaky waves for the structure listed in Table 2 are computed for λ = 1 and = 1/2, which 1008 is less than B . For Mode I, as illustrated in Fig. 5 (a), where γ 0 is the initial guess value of the propagation constant, 1048 and the matrix G γ is the derivative of the G matrix elements 1049 with respect to γ at γ = γ 0 . The next iteration point is where χ is the maximum eigenvalue of G −1 0 G γ . However, 1052 as the number of space harmonics increases, for a given node 1053 distribution, the matrix G 0 becomes ill conditioned because 1054 of the large imaginary parts of h ν in (34) and of h l n in (45). 1055 When h l n d l in the latter equation has a large imaginary part, 1056 the cosine and sine expressions are identical, for a given 1057 precision, and the transfer matrix is singular and the grand 1058 matrix (55) becomes ill-conditioned for γ values in regions 1059 near the root.

1060
The condition number of a matrix can be used as a gauge 1061 to determine if it is singular [34], [41], [42], [43]. When 1062 the condition number is approximately zero, the matrix is 1063 well conditioned or non singular, whereas if the condition 1064 number is large, the matrix is singular and its determinant 1065 is approximately zero. Contemporary linear algebra software 1066 is used here to estimate the reciprocal condition number, R c , 1067 of G. In particular, the reciprocal condition number of G will 1068 be driven from a value close to unity, for the initial guess of γ 0 , 1069 to a value close to zero, at the desired complex root.

1070
For the structure in Fig. 11 and a grating period = 0.62λ 1071 which is below the second Bragg condition, the initial value 1072 γ 0 /k = 0.0004 + j 1.57 yields a root γ /k = 0.0023668 + 1073 j 1.5812627 for Mode I. The roots for the other modes are 1074 −0.0023668 + j 1.5812627 for Mode II and ±0.00276343 + 1075 j 1.57210132 for Modes III and IV. Thus, the initial guess for 1076 γ 0 should produce a large value of R c for the grand matrix 1077 compared to the machine precision. Starting with a small 1078 number of space harmonics, the reciprocal condition number 1079 for the initial guess is illustrated in Fig. 12. 1080 The first iteration of the Newton search involves the cal-1081 culation of R c for the initial value γ 0 . For the uniform node 1082 spacing of 0.2 λ, and three space harmonics, the initial R c is 1083 about 10 −3 . Using Newton iterations (∼3-4), R c → 10 −20 1084 and the Newton step (for γ ) is less than 10 −12 which is close 1085 to machine precision. On the other hand, as the number of 1086 space harmonics is increased above ∼ 20 the initial computed 1087 value of R c drops below ∼ 10 −15 , and consecutive iterations 1088 are unsuccessful in reducing R c . Indeed, for the number of 1089 space harmonics above about 27, R c fluctuates and step sizes 1090 approach zero.

1091
When the node separation is 0.1 λ, initial values of R c are 1092 higher for identical numbers of space harmonics. Indeed, the 1093 lower node spacing produces stable results for a number of 1094 space harmonics exceeding 45. The cost for the additional 1095 accuracy is computation cycles.

1096
The process of adding space harmonics to obtain a Cauchy-1097 type convergence for γ is addressed below. The circle 1098  Table 1. Thus, to get proper root differences, 1119 two different cases were computed: (1) when the maximum 1120 space harmonic produces a zero value on the diagonal and

1125
This section will be concerned with properties of the waveg-1126 uide mode and its various space harmonics in the grating layer 1127 for the structure given in Table 2 and illustrated in Fig. 11. The  constant and the layer thickness. However, additional param-1135 eters such as the grating duty cycle affect the mode char-1136 acteristics. The material in the grating regions, specified in 1137 Table 2 represent the material that forms the tooth regions. 1138 The material of the fill regions is air.

1139
Calculations for Mode I are made for the structure given 1140 in Fig. 11 with seven nodes and a node spacing of 0.1 µm. 1141 The number of space harmonics are indicated in the various 1142 figures. Calculations for various modal characteristics are 1143 plotted in Figs. 14-17 as a function of (β − K )/K where the 1144 grating period will be varied and the complex propagation 1145 constant γ = α + jβ.   Computations for the three duty cycles show that Mode I, 1194 Fig. 18 (a), has a positive attenuation coefficient before the 1195 second Bragg condition but it changes sign after the second 1196 Bragg condition when > B . Mode III shows similar 1197 behavior for the 50 and 75% duty cycles but the 25% duty 1198 cycle case has α III < 0 for < B and α III > 0 for > B . 1199 In particular, when d c ≈ 0.3 , α III > 0 but α III ≈ 0 for < 1200 B . Fig. 19 shows Mode III attenuation for the duty cycles 25, 1201    into the superstrate. If the output power from the superstrate 1218 can be directed into the substrate, the mode will lossless. If the 1219 waveguide has a finite-depth substrate, then two waveguides 1220 could coupled by way of leaky mode radiation.

1221
The 30 and 29.9 percent duty cycle curves for Mode III 1222 of Fig. 21 also have attenuation coefficients that are zero at 1223 different grating periods < B . On the other hand, Mode I 1224 has α > 0 while Mode II has α < 0 at all values of < B 1225 provided β −1 lies in the fast-wave region.

1226
The sharp changes in α around the second Bragg condition 1227 (dα/d → ±∞) implies a very narrow bandwidth for either 1228 no radiation from the waveguide mode or for coupling power 1229 into the waveguide mode. of Modes II and IV would require the plane wave directed 1253 at the substrate, Mode IV, be at a different angle from that 1254 of Mode II with power headed toward the central waveguide 1255 region. When ≈ B , the substrate/superstrate angles of 1256 plane-wave radiation for Modes I and II become indistin-1257 guishable from those of Modes III and IV.

1258
For an infinite-length grating waveguide, all modes have 1259 α = 0 at the exact second Bragg condition because the 1260 amplitude of the n = −1 space harmonic tends to zero as 1261 → B . However, note that for the waveguide of Table 2, 1262 Fig. 9 indicates that the bandwidth about α = 0 is very small. 1263 Nevertheless, when = B , power can neither leak in nor 1264 leak out, i.e., power cannot be coupled in or out of the waveg-1265 uide mode. On the other hand, a section of a finite-length 1266 grating connecting longitudinally between two non-grating 1267 waveguides may provide coupling due to the excitation of 1268 radiation modes [44] induced at intersections of the a periodic 1269 and non-periodic waveguides.

1271
Practical applications use finite-length gratings. Periodic 1272 waveguide structures are used in various devices and systems 1273 such as reflectors for distributed feedback lasers, enhance-1274 ment of coupling between dielectric waveguides, and emis-1275 sion of radiation in a specific direction in the superstrate 1276 or substrate regions. Finite-length periodic waveguides are 1277 addressed below.

1278
The waveguide structure illustrated in Fig. 1 may be 1279 viewed as a collection of waveguides placed back-to-back. If 1280 materials in the tooth and fill regions along with all the other 1281 layers are lossless, then the modes in each of those regions 1282 consist of a set of discrete-spectra modes whose fields tend 1283 to zero as x → ±∞, resulting in a set of continuous-spectra 1284 modes (radiation modes) whose fields remain finite at all 1285 positions in the superstrate and substrate regions, along with 1286 a set of continuous-spectra evanescent modes confined to the 1287 vicinity of the discontinuity. These three sets of modes form 1288 a complete orthogonal set. Continuous-spectra modes are 1289 induced at the tooth and fill boundaries and thus scatter power 1290 into or out of the central waveguide region. The process of 1291 matching fields with discrete and continuous spectra, in the 1292 tooth region to the those in the fill region produces the fields 1293 in the periodic waveguide.

1294
Floquet theory, in essence, condenses the mode sets to 1295 one mode whose field has a periodicity and a propagation 1296 constant as given in (6). At grating periods near the second 1297 Bragg condition, there is a single space harmonic, n = −1, 1298 that has exponential growth in the claddings with propagation 1299 either in or out of the cladding regions. At points x 1 and 1300 x 5 in Layers 1 and 6 of Fig. 1, the Poynting vector indicates 1301 the direction of power flow due to the −1 space harmonic. 1302 Power computed from the Poynting vectors correlates with 1303 the power loss or gain of the computed value of α.

1304
In a typical application, a periodic waveguide is located 1305 longitudinally between waveguides without gratings. The 1306 fields of the input waveguide are characterized by incident 1307 are characterized only by transmitted fields. The fields in 1309 the periodic waveguide will have forward and backward Since the root search for the propagation constant is devel-1362 oped using Newton's method, one needs the derivative of the 1363 eigenvalues and eigenvectors with respect to γ , obtained from 1364 (29) and by the vector expansion given by (30). Furthermore, 1365 Appendix B describes a method of calculating the group 1366 velocity of the modes which involves the calculation of the 1367 derivative of the grand matrix with respect to γ and the deriva-1368 tive with respect to the free-space wavenumber k. An iteration 1369 scheme such as one that varies the grating period , after a 1370 root is determined, allows a new root at the iteration point 1371 to be approximated from the derivative dγ /d , obtained by 1372 differentiation of the grand matrix, (55), with respect to γ 1373 and . Note that there is a fourth-order zero at the second 1374 Bragg condition, so that the Newton root search has trouble. 1375 However, Figs. 9 and 17 illustrate sufficient accuracy in the 1376 neighborhood of the second Bragg condition. The expla-1377 nation of why second-order gratings are successfully used 1378 as reflectors and outcouplers in optical waveguide devices 1379 is explained physically based on the mathematical analysis 1380 developed in this paper.

1381
In addition to showing mathematically that there is no 1382 mode attenuation at the exact second Bragg condition, 1383 we show that some modes with leaky waves coming into the 1384 substrate and leaky waves going out of the superstrate have 1385 a second zero in the attenuation coefficient at grating periods 1386 near the second Bragg condition.

1387
There are two sources of dispersion in dielectric waveg-1388 uides: (1) material and (2) waveguide geometry. In most 1389 waveguide geometries the material dispersion plays a minor 1390 role if the operating wavelength is significantly different from 1391 the material's bandgap wavelength. Hence, the source of the 1392 derivative of the grand matrix, G k in Appendix III-B is due to 1393 the parameter k.

1395
An iterative method is used for driving the eigenvalue υ(γ , k) 1396 of the equation to zero, where the complex variables γ and k are the propaga-1399 tion constant and free-space wavenumber, respectively. For a 1400 given wavelength the matrix G and its associated eigenvec-1401 tor r are functions of γ . Similarly, for a given propagation 1402 constant, G and its eigenvector are functions of k.

1403
The object here is to determine the value of γ , for a 1404 specific value of k = k a that renders G as singular and 1405 thus maps the non-trivial eigenvector r to 0. If γ = γ a 1406 makes G singular, then r(γ a ) is the non-trivial solution vector, 1407 or eigenvector whose eigenvalue is 0. (There may be more 1408 than one solution.) If for an arbitrary value of γ , G is non-1409 singular, G cannot map r, whose norm is unity, to 0, but 1410 it can, for example, map the eigenvector associated with its 1411 VOLUME 10, 2022 eigenvector, so one possibility is to make the vector r the The group velocity can be calculated from the matrix G(γ , k) 1449 when it is singular. We outline the formulation below. That is, G H (γ a , k a ) maps the vector s a ( s a = 0) to 0.

1456
In this paper there are two types of eigen expansions: 1457 a) expansion for Newton's method for root searching, and 1458 b) expansion about singular points for evaluating wave prop-1459 erties by way of eigenvectors. For a) the expansion (A.2) 1460 is about the point γ = γ 0 and a fixed value of k, where 1461 G is non singular. For b) it is assumed that both k and γ 1462 are independent variables, and (A.1) will be expanded about 1463 k = k a and γ = γ a . When G, υ and r are expanded about the 1464 singularity, (A.1) becomes, keeping first two terms, 1465 G a + G k δk + G γ δγ r a + r k δk + r γ δγ 1466 = υ γ δγ + υ k δk r a + r k δk + r γ δγ , (B.2) 1467 where υ(γ a , k a ) = 0. The terms with subscripts of k and 1468 γ represent partial derivatives with respect to the variables 1469 evaluated at k = k a , and γ = γ a . The first-order terms in 1470 (B.2) become 1471 G a r k δk + G a r γ δγ + (G k δk + G γ δγ )r a where c is the velocity of light in vacuum. We defined mode 1488 propagation along the z direction as exp(jωt − γ z) where γ is 1489 complex. However, group velocity which specifies the veloc-1490 ity of a wave packet is usually defined as the real quantity 1491 1/(dβ/dω). In (B.7) we use the transformation γ = jβ so 1492 thatβ is complex.