1126 JOURNAL OF LIGHTWAVE TECHNOLOGY,VOL.6.NO.6.JUNE 1988 Guided Modes of Ti:LiNbO,Channel Waveguides: A Novel Quasi-Analytical Technique in Comparison with the Scalar Finite-Element Method E.STRAKE,G.P.BAVA.AND I.MONTROSSET Abstract-Two different techniques for the calculation of optical cially for Ti:LiNbO:channel waveguides.This very fast modes of Ti:LiNbO,channel waveguides are presented.The first one method is well suited for optimization procedures,where is an almost analytical technique based on the effective-refractive-in- dex method.The second one is based on the well-known finite-element numerous calculations may have to be performed in a rea- method.Both approaches use a rather realistic waveguide model con- sonable time.However,as it is based on the effective- sidering crystal anisotropy,wavelength dispersion of the refractive in- refractive-index method (ERI)and on suitable approxi- dex distributions,and also the nonlinear dependence of the ordinary mations,some errors are introduced having a magnitude refractive index on the Ti'concentration.Effective indices and field hard to estimate.To assess the quality of the new tech- distributions obtained by both techniques are compared for wave- nique and to quantify the induced errors,we further use lengths of 0.6 and 1.2 um,for both quasi-TE and-TM polarization, considering Y-cut and Z-cut crystal orientation and both isotropic and the purely numerical finite-element method(FEM)which anisotropic Ti diffusion. yields highly accurate results.Then a detailed comparison allows us to recognize the limits of validity of the quasi- I.INTRODUCTION analytical approach. The two methods are used with a rather realistic model PTICAL waveguides in LiNbO3 crystals,obtained by means of titanium in-diffusion.form the basis of of the waveguiding structure and applied to both Z-cut and Y-cut (or X-cut)configurations by including the effects of many integrated optical devices.Due to several physical properties of the crystal material,e.g.,low loss,large crystal anisotropy,nonlinear dependence of the ordinary electrooptical,piezoelectrical,and elastooptical coeffi- refractive index on the local Ti concentration.aniso- cients and high second-order nonlinearity,a large number tropic Tidiffusion and wavelength dispersion of the re- of interesting devices can be realized. fractive index distributions.The waveguide model is pre- The optimization of the efficiency of integrated optic sented in Section II. devices requires the knowledge of the effective refractive Sections llI and IV describe the two computational indices and field distributions of the waveguiding struc- methods in detail. ture and their dependence on the fabrication parameters Examples of the results are given in Section V for two Therefore,a great interest in theoretical methods of wave- wavelengths (0.6 and 1.2 um).both polarizations (quasi- guide analysis exists. TE and quasi-TM),two crystal orientations(Y-cut and Z- Several papers,using different approaches and different cut),and both isotropic and anisotropic Ti diffusion. degrees of approximations,have been concerned with the Field distributions and effective refractive indices ob- electromagnetic analysis of anisotropic and diffused chan- tained by the two methods are compared in detail and pos- nel waveguides (e.g..[1]-[4]).However,information sible fields of application are discussed about the comparison of the results obtained by different techniques,especially for the field distributions,can II.WAVEGUIDE MODEL hardly be found. A.Diffusion Profiles It is the aim of this paper to present first a novel quasi- analytical technique of solving the wave equation espe- The assumed configuration with a representation of the waveguide axes is shown in Fig.1.Optical channel wave- Manuscript received July 17.1987:revised November 6.1987. guides can be formed in LiNbO;by the in-diffusion of E.Strake is with Angewandte Physik.Universitat-GH-Paderborn.4790 thin titanium stripes.On the basis of a widely used dif- Paderborn.W.Germany. fusion model and taking into account anisotropy effects, G.P.Bava is with Dipartimento di Elettronica,Politecnico di Torino. 10129 Torino.Italy. the local Ti concentration c(x,y)can be written as 151. 1.Montrosset is with DIBE.Universita di Genova.16145 Genova.It-6]: aly IEEE Log Number 8718899. c(x.y)=cof(u)g(s) (1) 0733-8724/88/0600-1126$01.00©1988IEEE
1126 JOURNAL OF IIGHTWAVF TFCHNOLOGY VOL 6 NO 6 JUNF 1Y88 Guided Modes of Ti : Limo3 Channel Waveguides: A Novel Quasi-Analytical Technique in Comparison with the Scalar Finite-Element Method E. STRAKE. G. P. BAVA. AND I. MONTROSSET Abstract-Two different techniques for the calculation of optical modes of Ti : LiNb03 channel waveguides are presented. The first one is an almost analytical technique based on the effective-refractive-index method. The second one is based on the well-known finite-element method. Both approaches use a rather realistic waveguide model considering crystal anisotropy, wavelength dispersion of the refractive index distributions, and also the nonlinear dependence of the ordinary refractive index on the Ti'+ concentration. Effective indices and field distributions obtained by both techniques are compared for wavelengths of 0.6 and 1.2 pm, for both quasi-TE and -TM polarization, considering Y-cut and Z-cut crystal orientation and both isotropic and anisotropic Ti" diffusion. I. INTRODUCTION PTICAL waveguides in LiNb03 crystals, obtained by 0 means of titanium in-diffusion, form the basis of many integrated optical devices. Due to several physical properties of the crystal material, e.g., low loss, large electrooptical, piezoelectrical, and elastooptical coefficients and high second-order nonlinearity, a large number of interesting devices can be realized. The optimization of the efficiency of integrated optic devices requires the knowledge of the effective refractive indices and field distributions of the waveguiding structure and their dependence on the fabrication parameters. Therefore, a great interest in theoretical methods of waveguide analysis exists. Several papers, using different approaches and different degrees of approximations, have been concerned with the electromagnetic analysis of anisotropic and diffused channel waveguides (e.g., [ I]-[4]). However, information about the comparison of the results obtained by different techniques, especially for the field distributions, can hardly be found. It is the aim of this paper to present first a novel quasianalytical technique of solving the wave equation espeManuscript received July 17, 1987: revised November 6, 1987. E. Strake is with Angewandte Physik, Universit~t-GH-Paderborn. 4790 G. P. Bava is with Dipartimento di Elettronica, Politecnico di Torino. I. Montrosaet is with DIBE. Universita di Genova, 16145 Genova. ItIEEE Log Number 8718899. Paderborn. W. Germany. 10129 Torino. Italy. aly. cially for Ti : LiNb03 channel waveguides. This very fast method is well suited for optimization procedures, where numerous calculations may have to be performed in a reasonable time. However, as it is based on the effectiverefractive-index method (ERI) and on suitable approximations, some errors are introduced having a magnitude hard to estimate. To assess the quality of the new technique and to quantify the induced errors, we further use the purely numerical finite-element method (FEM) which yields highly accurate results. Then a detailed comparison allows us to recognize the limits of validity of the quasianalytical approach. The two methods are used with a rather realistic model of the waveguiding structure and applied to both Z-cut and Y-cut (or X-cut) configurations by including the effects of crystal anisotropy, nonlinear dependence of the ordinary refractive index on the local Ti4+ concentration, anisotropic Ti4+ diffusion and wavelength dispersion of the refractive index distributions. The waveguide model is presented in Section 11. Sections I11 and IV describe the two computational methods in detail. Examples of the results are given in Section V for two wavelengths (0.6 and 1.2 pm), both polarizations (quasiTE and quasi-TM), two crystal orientations ( Y-cut and Zcut), and both isotropic and anisotropic Ti4+ diffusion. Field distributions and effective refractive indices obtained by the two methods are compared in detail and possible fields of application are discussed. 11. WAVEGUIDE MODEI A. Diffusion Projiles The assumed configuration with a representation of the waveguide axes is shown in Fig. 1. Optical channel waveguides can be formed in LiNbOl by the in-diffusion of thin titanium stripes. On the basis of a widely used diffusion model and taking into account anisotropy effects, the local Ti4+ concentration c(x, y) can be written as [5], [GI: c(x, .VI = cof(u) g(s) (1) 0733-8724/88/0600-1126$01 .OO O 1988 IEEE
STRAKE et al.:GUIDED MODES OF Ti:LiNbO,CHANNEL WAVEGUIDES 1127 TLiNbO3 LNb03 Fig.1.Ti:LiNbO,waveguide configuration. 00.20406Q8 10 Co[102 cm3] with Fig.2.Ordinary and extraordinary refractive index dependence on the Ti concentration,Comparison between the experimental(●，×)and the f(u)=exp(-u), analytical dependences used in our model. 8(s) expressions were adopted: erf ((W/2D,)(1 s))erf ((W/2D,)(1 -s)) 0.839(入/um)月 d(入)= 2 (入/um）-0.0645 (2) d,(入)= 0.67(入/um)月 (6) where u =y/D,and s 2x/W are normalized coordi- (入/um)-0.13 nates,D.and D.are the diffusion depths along the two A comparison between measured values [9]and the transverse directions,W is the Ti stripe width before dif- previous analytical expressions for on.(c)is presented fusion,and co is a parameter connected to the total Ti-ion in Fig.2.A small discrepancy can be noted for on at low content per unit length along z,given by coD,W/2.values of ci for on,using F=1.3 x 10-25 cm'and y= By introducing the atomic weight and bulk density of 0.55,discrepancies arise forc>0.9x102 cm-3.The Ti.we can relate co to the Ti stripe thickness 7 before waveguides considered in this paper belong to co=0.5 diffusion [7]: x 1021 cm3 where the difference between on,and on is T=acoD. (3) well approximated by (5). with a=1.57×1023cm3. B.Wave Equation The variation of the bulk crystal refractive indices with Assuming a properly oriented crystal exhibiting a di- temperature and wavelength is considered using Sellmeier agonal permittivity tensor,introducing the usual approx- equations given in [8]. imation of quasi-TE and quasi-TM modes [13],and ne- The local increase of the extraordinary (on(x,y))and glecting the transverse refractive index gradient as ordinary (on (x,y))refractive indices with respect to the compared with the field gradient,a scalar,two-dimen- bulk values (n and n respectively)is related to c by sional wave equation for the dominant field component can be obtained [14]: 6n,(入.c）=d(入)f(c),i=o,e (4) where f gives the dependence of the refractive index vari- +a,2 n+knxy)）-g:》 (x,y)=0 ation on c at a given wavelength and d,(A)takes into account the effects of dispersion.Using the results in [9], (7) [10]at A =0.6328 um as reference,we obtained the fol- where ko =2/A is the free-space wavenumber.B lowing approximating expressions: koNerr is the mode propagation constant,and (x.y)is f(c)=Ec the transverse field distribution (E.for quasi-TE and H for quasi-TM modes )The relations between the wave f(c)=(Fc)". (5) equation parameters n and the field for the two polarizations and the elements n.nn characterizing the In (5)a linear relation is assumed between on,and c(E permittivity tensor are reported in Table 1.In Table Il the =1.2 x 10-23 cm3),while a good fitting of the experi- relations between n,n.n:.D..D..and n.n.D.D mental on versus c dependence requires a nonlinear func- (diffusion depths along and perpendicular to the optical tion f(c).The choice of F and y is such to achieve good axis)for the three possible c-axis orientations are pre- agreement with the experimental data in the c-range of sented. interest.The dispersion of the refractive index changes was approximated by means of a simple oscillator model III.APPROXIMATE ANALYTICAL SOLUTION [9]in the infrared region:after a comparison of the results The approximate analytical technique is based on the available in the literature [9],[11],[12]the following well-known effective-refractive-index method [13].[15]
STRAKE ('I ul GUIDED MODES OF TI LINbO, CHANNEL WAVEGUIDES 1127 Fig. 1. Ti : LiNbOI waveguide configuration with f(u) = exp (-U?), g(s) - erf((W/2D1)(1 + 4) + erf((W2D.m - 4) - 2 (2) where U = y/D, and s = 2x/W are normalized coordinates, D, and D,. are the diffusion depths along the two transverse directions, W is the Ti stripe width before diffusion, and eo is a parameter connected to the total Ti-ion content per unit length along z, given by &coD,. W/2. By introducing the atomic weight and bulk density of Ti. we can relate co to the Ti stripe thickness T before diffusion [7]: T = UC~D, (3) with a = 1.57 x lopz3 cm3. The variation of the bulk crystal refractive indices with temperature and wavelength is considered using Sellmeier equations given in [8]. The local increase of the extraordinary ( 6nl,(x, y)) and ordinary ( 6n, (x, y ) ) refractive indices with respect to the bulk values ( nch and nob, respectively) is related to c by 6ni(X, c) = dj(X)i(c), i = 0, e (4) where j; gives the dependence of the refractive index variation on c at a given wavelength and dj (A) takes into account the effects of dispersion. Using the results in [9], [IO] at X = 0.6328 pm as reference, we obtained the following approximating expressions: fc(c) = Ec fo(c> = (Fc)'. (5) In (5) a linear relation is assumed between 6n,, and c( E = 1.2 X cm3), while a good fitting of the experimental 6n, versus c dependence requires a nonlinear functionL,( e). The choice of F and y is such to achieve good agreement with the experimental data in the c-range of interest. The dispersion of the refractive index changes was approximated by means of a simple oscillator model [9] in the infrared region; after a comparison of the results available in the literature [9], [ll], [12] the following m 05 c" 3 0 !:E 0 0.2 0.L 0.6 0.8 1.0 ~~110~~~~1 - Fig. 2. Ordinary and extraordinary retractive index dependence on the Ti concentration. Comparison between the experimental ( 0. X ) and the analytical dependences used in our model. expressions were adopted: 0.839 (A/pm>' (h/prn)? - 0.0645 0.67 (Xipm)' (Xjprn)' - 0.13' d,,(h) = d,,(X) = A comparison between measured values [9] and the previous analytical expressions for 6n,,,,, (c) is presented in Fig. 2. A small discrepancy can be noted for 6n,, at low values of c;-for 6n,, using F = 1.3 x IO-'' cm' and y = 0.55, discrepancies arise for c > 0.9 x IO" cmP3. The wave uides considered in this paper belong to c0 = 0.5 X 10" cmp3 where the difference between 6n,, and 6n,, is well approximated by (5). B. Wave Equation Assuming a properly oriented crystal exhibiting a diagonal permittivity tensor, introducing the usual approximation of quasi-TE and quasi-TM modes [ 131, and neglecting the transverse refractive index gradient as compared with the field gradient, a scalar, two-dimensional wave equation for the dominant field component can be obtained [ 141: (7) where ko = 21r/X is the free-space wavenumber, 0 = koN,,, is the mode propagation constant, and $(x, y) is the transverse field distribution (E, for quasi-TE and H, for quasi-TM modes). The relations between the wave equation parameters a,, a\, n and the field $ for the two polarizations and the elements n,, n,, n; characterizing the permittivity tensor are reported in Table 1. In Table I1 the relations between n,, n,, n;, D,, D,, and n,,, n,,, D,,, D, (diffusion depths along and perpendicular to the optical axis) for the three possible c-axis orientations are presented. 111. APPROXIMATE ANALYTICAL SOLUTION The approximate analytical tcchnique is based on the well-known effective-refractive-index method [ 131, [ 151
1128 JOURNAL OF LIGHTWAVE TECHNOLOGY.VOL..6.NO.6.JUNE 1988 TABLE I Optica】field polarization quasi-TE 2 1 quasi-TM 02y TABLE II Relation between Usual configuration crystal coordinate nanes nx Dx X-cut,Y-propagation c-axis× 9 Y-cut,X-propagation Z-cut,X-propagation c-axis ll y no ne Z-cut,Y-propagation X-cut,Z-propagation c-axis I z e 9 Y-cut,Z-propagation and assumes the parameters a and o,to be constant in the waveguide cross section and evaluated by inserting 10 the bulk refractive index values (cf.,(7)and Table I).So the waveguide is completely characterized by n(x,y). A.Field Computation As a first step,we separate(7)into two equations.By n88 putting (s,u)=X(s)Y(u;s),we get approximately a'y 05 (koD)((s.u)())Y=0 (a) a80 d'x +(a,kw/2y'(a(s)-Nm)x=0.(8b) The semicolon in Y(u:s)points out that s plays the role of a parameter in (8a),while (s)is a local eigen- 0 value which,when used in(8b),allows to obtain the mode 05101520 effective index Nm. In (8a)one has n'(s,u)=ng+2n6n(s.u),where np Fig.3.Comparison between the ordinary and extraordinary refractive in- is the bulk refractive index value.In order to proceed to em7osg60a"0o88 10cm3and入=0.6328 an analytical solution of (8a)we approximate on(s,u)as follows [16],[17]: on(s.u)=on(s,0)(1-M tanh2 (nu))(9) By introducing p(s)=2n 6n(s.0)and changing from where M and n are two parameters.Some examples using the variable u tog=inu+/2.(8a)can be transformed this approximation are shown in Fig.3 for the case M into: I and M=1.05 for the ordinary and extraordinary in- a'Y dices.The corresponding values of n are 0.820 and 0.776 02 +(a,kD,/n）（nm- for the ordinary,1.085 and 1.035 for the extraordinary index.In the figure on is normalized with respect to the on (0,0)value at A=0.6328 um. +(M-1)p- Mp Y=0 (10) sin
I128 JOURNAL OF LIGHTWAVE TECHNOLOGY. VOI. 6. NO 6. JUNE IYXX quasi-TE TABLE I Optical field n - ;z 1 EX X Relation between I/ n D crystal coordinate system and waveguide nx coordinate system Y nZ Dx Y c-axis 11 x n n n D/~ Dl c-axis 11 y n n "o D_L '11 "e D_L D_L c-axis 11 z Usual configuration names X-cut, Y-propagation Y-cut, X-propagation Z-cut. X-propagation Z-cut, Y-propagation X-cut, 2-propagation Y-cut, Z-propagation and assumes the parameters a, and a, to be constant in the waveguide cross section and evaluated by inserting the bulk refractive index values (cf., (7) and Table I). So the waveguide is completely characterized by n(x, y). A. Field Computation putting $(s, U) = X(s) Y(u; s), we get approximately As a first step, we separate (7) into two equations. By d'X ~ + (c~,k~W/2)~ (&(s) - N&) X = 0. ds (8b) The semicolon in Y(u; s) points out that s plays the role of a parameter in (8a), while &(s) is a local eigenvalue which, when used in (8b), allows to obtain the mode effective index Ne,,. In (8a) one has n2(s, U) = ni + 2nh6n(s, U), where nh is the bulk refractive index value. In order to proceed to an analytical solution of (8a) we approximate 6n (s, U) as follows [16], [17]: 6n(.s, U) = 6n(s, 0)(1 - Mtanh2 (vu)) (9) where M and 7 are two parameters. Some examples using this approximation are shown in Fig. 3 for the case M = 1 and M = 1.05 for the ordinary and extraordinary indices. The corresponding values of 7 are 0.820 and 0.776 for the ordinary, 1.085 and 1.035 for the extraordinary index. In the figure 6n is normalized with respect to the 6n,(0, 0) value at X = 0.6328 pm. 05 10 15 20 u- Fig. 3. Comparison between the ordinary and extraordinar) refractive index distributions (-) in depth and our approximation with M = I .O (. . . ) and M = 1.05 (---) for c,, = 0.3 . IO" cin' and X = 0.6328 w. By introducing p (s) = 2 tzI, 6n (s, 0) and changing from the variable U to t; = i7u + ~/2, (8a) can be transformed into: at;
STRAKE e 1.:GUIDED MODES OF Ti:LiNbO,CHANNEL WAVEGUIDES 1129 which is of the form striction because the ERI-method is in principle limited to well guided modes,due to the approximate nature of +(g+4')2+'(1-) o'y Y=0.(11) the decoupling of the two-dimensional wave equation (7) sin into (8a)and (8b). Restricting the analysis to guided modes and assuming B.Normalization the approximate boundary condition Y=0 for u=0,we obtain solutions of (11)for odd values of g'(g'=2g Using the previous results,the expression for the trans- I,q=0,l,2.···),that can be expressed as[18: verse electric or magnetic field component (E,for TE,H. for TM)is: Y()=(sin)“C+,(cos) (12) (s,u)=Npoxp(s)Y(u:s) where C are Gegenbauer polynomials as functions of cos). =Npg(cosh ks)"Cr(-i sinh Ks) The assumed boundary condition is a rather good ap- proximation owing to the large refractive index difference (cosh nu)"C(-i sinh nu).(17) between the LiNbO,crystal and the cover material (in Its normalization according to the relation [19] most cases the cover material is air).With this assumption Y,(;s）=(cosh nu)Ci(-i sinh nu）(13) 2P= (x)2.d and,by comparison of(10)and (11),we obtain: leads to the following approximate expression of N u'(s)=(1-V1+4Mp(s)(a,kD,/n)) D.WNer 4PT Jx(s nr(s)=n呢+(1-M)p(s) +(n/akD)(q'+u'(s)) 14) (u:s)du ds (18 In (13).the number g defines the mode order in depth and Yo(u;s)is the depth distribution of the fundamental where T Zo for TE-modes and T =1/Zo for TM-modes: mode.The'effective index profile''ne(s)given by (14) Zo is the free space impedance and P denotes the optical must be used successively in the solution of the lateral power carried by the mode. field equation (8b). Equation (18)can be evaluated analytically by consid- In order to find X(s)from the solution of(8b),a similar ering Y(u:s)=Y(u;0).Then the integral factorizes approximate procedure is applied by fitting ner(s)with and one obtains,e.g.,for the fundamental TE mode [17], the function [181: nr(s）=i.+δnem(1-R tanh2(ks)）(15) INoo2=D.WNem -2-2a+a1-5r2(-μ) where onm=n()-n and with a suitable choice of nKZo T(-2) the parameters R and K.A treatment completely similar to the case of (8a)leads to: T---4-2 (19) T(-2μ'-2)r(-2μ')/ X(s)=(cosh Ks)"Co(-i sinh Ks) where u'='(0)is obtained from (14)and I is the =(11+4R(a,ko Woncto/2K)) gamma function [18].The result(19)is quite useful since it allows one to avoid the numerical integration needed Ne ni onc (1-R)+(2K/askoW)(p +u) for field normalization. (16） IV.FINITE-ELEMENT SOLUTION where p=O,l,2.···gives the mode order in the lateral One of the methods for solving the two-dimensional direction wave equation without the need of additional approxi- For guided modes the field distribution vanishes for mations is the finite-element method,as described,for in- large values of s or u.From (13)and (16)it follows stance,by [2].[3].[14],[20],and [21].This method al- that (p+u)and (q'+')must be negative.If one of lows to solve the wave equation exactly also in the these numbers becomes zero,the mode undergoes'cut- vectorial form.Several authors [3],[22]compared dis- off.Usually this is expressed in terms of the effective persion curves for rectangular dielectric waveguides,ob- index as Ner lcutor =n tained by scalar and vectorial FEM calculations and found The case M>I or R>1,which may occur in order a very good agreement.Ti:LiNbO3 waveguides having to obtain a good index profile approximation near the very small refractive index changes without lateral dis- waveguide axis,leads to cutoff effective index values de- continuities are much more similar to planar waveguides viating from n.This means,however,no additional re- for which the scalar formulation is rigorous.The FEM
STRAKE [’f ol GUIDbD MODES OF TI LiNbO, CHANNEL WAVEGUIDES 1129 which is of the form striction because the ERI-method is in principle limited to well guided modes, due to the approximate nature of the decoupling of the two-dimensional wave equation (7) into (8a) and (8b). Y = 0. ( 11 Restricting the analysis to guided modes and assuming the approximate boundary condition Y = 0 for u = 0, we obtain solutions of (1 1) for odd values of q’( qr = 2 q + 1, q = 0, 1, 2, * * e), that can be expressed as [lS]: ~(4) = (sin 4)’’ cos c;) (12) where C:/v’) are Gegenbauer polynomials (as functions of cos 4). The assumed boundary condition is a rather good approximation owing to the large refractive index difference between the LiNbO, crystal and the cover material (in most cases the cover material is air). With this assumption Y[,(u; s) = (cosh qu)” Cig’l I (-i sinh vu) (13) and, by comparison of (10) and (1 l), we obtain: B. Normalization Using the previous results, the expression for the transverse electric or magnetic field component (E, for TE, H, for TM) is: 4(s, = N/,yxp(s) Yyb; s) = N,,,(cosh KS)’ CjIp)( -i sinh KS) (cosh vu)” CiG’i I ( -i sinh vu). Its normalization according to the relation [ 191 (17) +m 2P = ST_ 1 (E x li*) . e’;dxdy --m leads to the following approximate expression of In (1 3), the number q defines the mode order in depth and Y,,(u; s) is the depth distribution of the fundamental mode. The “effective index profile” nzrf( s) given by (14) must be used successively in the solution of the lateral field equation (Sb). In order to find X( s) from the solution of (Sb), a similar approximate procedure is applied by fitting n,Zff( s) with the function n$,(s) = 12: + 6rz&,( 1 - R tanh’ (KS)) (15) where 6n& = &(O) - ni and with a suitable choice of the parameters R and K. A treatment completely similar to the case of (Sa) leads to: X,,(s) = (cosh KS)’ CL”( -i sinh KS) p = ~(1 2 - JI + 4~(a,k~~6n,,,/2K)~) 2 N:,, = n; + 6n:fr,,(1 - R) + (2K/a.,kow)’ (p + p) (16) wherep = 0, I, 2, . . * gives the mode order in the lateral direction. For guided modes the field distribution vanishes for large values of 1 s 1 or I U 1. From (13) and (16) it follows that (p + p) and (q’ + p’) must be negative. If one of these numbers becomes zero, the mode undergoes “cutoff’. Usually this is expressed in terms of the effective index as N& (cutoff = ni. The case M > 1 or R > 1, which may occur in order to obtain a good index profile approximation near the waveguide axis, leads to cutoff effective index values deviating from n/,. This means, however, no additional rewhere T = Zo for TE-modes and T = 1 /Zo for TM-modes; Zo is the free space impedance and P denotes the optical power carried by the mode. Equation (1 8) can be evaluated analytically by considering Y,(u; s) = Y,(u; 0). Then the integral factorizes and one obtains, e.g., for the fundamental TE mode [ 171, where pr = p’(0) is obtained from (14) and r is the gamma function [lS]. The result (19) is quite useful since it allows one to avoid the numerical integration needed for field normalization. IV. FINITE-ELEMENT SOLUTION One of the methods for solving the two-dimensional wave equation without the need of additional approximations is the finite-element method, as described, for instance, by [2], [3], [14], [20], and [21]. This method allows to solve the wave equation exactly also in the vectorial form. Several authors (31, [22] compared dispersion curves for rectangular dielectric waveguides, obtained by scalar and vectorial FEM calculations and found a very good agreement. Ti : LiNb03 waveguides having very small refractive index changes without lateral discontinuities are much more similar to planar waveguides, for which the scalar formulation is rigorous. The FEM
1130 JOURNAL OF LIGHTWAVE TECHNOLOGY.VOL.6.NO.6.JUNE 1988 results can therefore be considered the reference solutions interpolated continuously.This can be achieved by intro- to the present waveguide problem. ducing the''nodal shape functions''[23]S.with the fol- Starting from (7)which we wish to solve for a wave- lowing properties: guiding structure described by (1)-(6).we proceed in three steps.First,the differential formulation is transformed Sp(xa yq)=8pg into an integral formula.Then the waveguide cross sec- tion is divided into discrete area elements with a poly- 0≤Sn(x,y)≤1 nomial approximation of the field distribution inside each S0,only within elements containing element.This yields an algebraic equation,which can then node p. be solved by well-known numerical methods. (25) A.Variational Formulation As shown,e.g..in [23],the following variational 6p is the Kronecker delta symbol and (x y)are the co- expression becomes stationary if a true solution of the ordinates of node g.This means that S is localized around wave equation (7)is inserted for the trial function node p and vanishes outside a certainneighborhood'of node p,depending on the chosen mesh division. Thus the field distribution is approximated by the se- ries: +k(Nm-n2)2）ddy (20 (x,y)= ∑nS,(x,y) (26) p= where with expansion coefficients being identical with the no- w=1, for quasi TE-modes dal field values: w=n2. for quasi-TM modes (21) pn=(xp,p） (27) and n,are defined in Table I. If we now insert (26)into (20).we obtain The characteristic equation.which is equivalent to the wave equation(7),is simply: L()= 名，三（-aw+N8bw), (28) g=l 6L=0 (22) with where 6''means the variation with respect to o. At the waveguide-cover interface (y 0)where the kn2S,5,- 1 asp as。. refractive indices are discontinuous,the correct behavior a:ax ax ofφis given by 1 aSp as continuous (23) (29) a dy dy dx dy wa中 continuous. (24) and a dy It is a great advantage of the variational formulation wkSpS dx dy (30) that due to the so-called natural boundary conditions [23], equation (24)is fulfilled automatically if is forced to Here means the cross section of element "'e'',and satisfy (23),which can be maintained most easily by a the sums extend over all elements proper choice of the trial function.Therefore even piece- Now (22)turns into a set of m equations wise linear trial functions are allowed without introducing additional complications. aL =0,p=1,2,··,m (31) B.Discretization Procedure don Next the x-y plane is divided into a number of subre- and can be written(with (28))in the usual matrix notation gions,the so-called finite elements.Within each of them the trial function is approximated by a suitably chosen A·b=NB·西 (32) polynomial.In the simplest case the elements are trian- gular and first degree polynomials are used.By this pro-where A and B are real symmetric matrices of order m. cedure the x-y plane is covered with a grid of discrete with elements given by ape and bp respectively.is a nodes,which coincide with the corners of the triangles. column vector containing the discrete field valuespp The continuous function is replaced by a set of dis- =1,2,···,m}. crete values{φp,p=l,2,·,,m},where m is the Equation(32)is an algebraic eigenvalue equation with total number of nodes.Between the discrete nodes,is eigenvalues Ner and eigenvectors o which correspond to
I130 JOURNAL OF LIGHTWAVE TECHNOLOGY. VOL. 6. NO. 6. JUNE 1988 results can therefore be considered the reference solutions to the present waveguide problem. Starting from (7) which we wish to solve for a waveguiding structure described by (1)-(6), we proceed in three steps. First, the differential formulation is transformed into an integral formula. Then the waveguide cross section is divided into discrete area elements with a polynomial approximation of the field distribution inside each element. This yields an algebraic equation, which can then be solved by well-known numerical methods. A. Variational Formulation As shown, e.g., in [23], the following variational expression becomes stationary if a true solution of the wave equation (7) is inserted for the trial function 4: + ki(N:f, - n’) 02) dx dy where w = 1, for quasi TE-modes , for quasi-TM modes (21) w 1 n-2 and a,,, a,., n, 4 are defined in Table I. wave equation (7), is simply: The characteristic equation, which is equivalent to the 6L = 0 (22) where “6” means the variation with respect to 6. At the waveguide-cover interface (y = 0) where the refractive indices are discontinuous, the correct behavior of 4 is given by 4, continuous (23 1 _- ” continuous. a: ay’ It is a great advantage of the variational formulation that due to the so-called natural boundary conditions [23], equation (24) is fulfilled automatically if 6 is forced to satisfy (23), which can be maintained most easily by a proper choice of the trial function. Therefore even piecewise linear trial functions are allowed without introducing additional complications. B. Discretization Procedure Next the x-y plane is divided into a number of subregions, the so-called finite elements. Within each of them the trial function is approximated by a suitably chosen polynomial. In the simplest case the elements are triangular and first degree polynomials ate used. By this procedure the x-y plane is covered with a grid of discrete nodes, which coincide with the corners of the triangles. The continuous function 4 is replaced by a set of discrete values { ’,,, p = 1, 2, . . , m}, where m is the total number of nodes. Between the discrete nodes, 4 is interpolated continuously. This can be achieved by introducing the “nodal shape functions” [23] S,, with the following properties: SP(XS7 Y,> = 4Iq 0 I S,,(X, y) I 1 S, f 0, only within elements containing node p. (25) S, is the Kronecker delta symbol and (x,, y,) are the coordinates of node q. This means that s,, is localized around node p and vanishes outside a certain “neighborhood” of node p, depending on the chosen mesh division. Thus the field distribution is approximated by the series: m ‘(x, y) = c ‘/,S,,(X? v) (26) p= I with expansion coefficients ‘/J being identical with the nodal field values: ’,, = ‘(xp Y/J. (27) (28) If we now insert (26) into (20), we obtain m 111 L(’) = c c (-apq + N,Z,b/J,) ‘p’q p= 1 ,= I with and PP Here !de means the cross section of element “e”, and Now (22) turns into a set of m equations the sums extend over all elements. and can be written (with (28)) in the usual matrix notation where A and B are real symmetric matrices of ord5r m, with elements given by aPy and bpq, respectively. 4 is a column vector containing the discrete field values { ‘JJ, p = 1, 2, , m}. Equation (32) is an algebraic eigznvalue equation with eigenvalues Nfff and eigenvectors C#I which correspond to