《光波导理论与技术 Optical Waveguides Principles and Technologies》课程教学资源（参考文献）Guided Modes of Ti：LiNbO3 Channel Waveguides：A Novel Quasi-Analytical Technique in Comparison with the Scalar Finite-Element Method

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-in￾dex method. The second one is based on the well-known finite-element method. Both approaches use a rather realistic waveguide model con￾sidering crystal anisotropy, wavelength dispersion of the refractive in￾dex 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 wave￾lengths 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 coeffi￾cients 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 struc￾ture and their dependence on the fabrication parameters. Therefore, a great interest in theoretical methods of wave￾guide analysis exists. Several papers, using different approaches and different degrees of approximations, have been concerned with the electromagnetic analysis of anisotropic and diffused chan￾nel 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 quasi￾analytical technique of solving the wave equation espe￾Manuscript 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. It￾IEEE 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 rea￾sonable time. However, as it is based on the effective￾refractive-index method (ERI) and on suitable approxi￾mations, some errors are introduced having a magnitude hard to estimate. To assess the quality of the new tech￾nique 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 quasi￾analytical 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, aniso￾tropic Ti4+ diffusion and wavelength dispersion of the re￾fractive index distributions. The waveguide model is pre￾sented 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 (quasi￾TE and quasi-TM), two crystal orientations ( Y-cut and Z￾cut), and both isotropic and anisotropic Ti4+ diffusion. Field distributions and effective refractive indices ob￾tained by the two methods are compared in detail and pos￾sible 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 wave￾guides can be formed in LiNbOl by the in-diffusion of thin titanium stripes. On the basis of a widely used dif￾fusion 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 coordi￾nates, D, and D,. are the diffusion depths along the two transverse directions, W is the Ti stripe width before dif￾fusion, 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 vari￾ation 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 fol￾lowing 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 experi￾mental 6n, versus c dependence requires a nonlinear func￾tionL,( 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 di￾agonal permittivity tensor, introducing the usual approx￾imation of quasi-TE and quasi-TM modes [ 131, and ne￾glecting the transverse refractive index gradient as compared with the field gradient, a scalar, two-dimen￾sional 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 pre￾sented. 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 eigen￾value 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 in￾dices. 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 in￾dex 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 ap￾proximation 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 trans￾verse 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 “cut￾off’. 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 de￾viating from n/,. This means, however, no additional re￾where 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 consid￾ering 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 approxi￾mations is the finite-element method, as described, for in￾stance, by [2], [3], [14], [20], and [21]. This method al￾lows to solve the wave equation exactly also in the vectorial form. Several authors (31, [22] compared dis￾persion curves for rectangular dielectric waveguides, ob￾tained by scalar and vectorial FEM calculations and found a very good agreement. Ti : LiNb03 waveguides having very small refractive index changes without lateral dis￾continuities 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 wave￾guiding structure described by (1)-(6), we proceed in three steps. First, the differential formulation is transformed into an integral formula. Then the waveguide cross sec￾tion is divided into discrete area elements with a poly￾nomial 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 piece￾wise linear trial functions are allowed without introducing additional complications. B. Discretization Procedure Next the x-y plane is divided into a number of subre￾gions, 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 trian￾gular and first degree polynomials ate used. By this pro￾cedure 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 dis￾crete 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 intro￾ducing the “nodal shape functions” [23] S,, with the fol￾lowing 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 co￾ordinates 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 se￾ries: m ‘(x, y) = c ‘/,S,,(X? v) (26) p= I with expansion coefficients ‘/J being identical with the no￾dal 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

STRAKE er al.:GUIDED MODES OF Ti:LiNbO,CHANNEL WAVEGUIDES 1131 the propagation constant and field distributions of optical modes in the waveguide channel. An important fact is that A and B depend on the wave- length A,but not on the effective refractive index Nerr. Thus the effects of wavelength dispersion (cf.,equation (6))can be included very easily.This is an additional ad- vantage of the scalar formulation as compared to vectorial techniques where,in general,the matrices depend on the effective refractive index while the wavelength is given by the eigenvalue.Then,in the vectorial formulation,due Fig.4.Example of a typical grid used for the FEM computations. to wavelength dispersion of the refractive indices of the waveguide,the matrices will also depend on A.so Nerr and The coordinate-overrelaxation method is an iterative A have to be specified as input parameters.As a result the technique of solving (32);therefore,an initial guess for eigenvalue will yield a wavelength value which,in gen- the field vector o is required.The number of iterations eral,will be different from the input one.In order to ob- needed to achieve convergence strongly depends on the tain a self-consistent vectorial solution the calculation has choice of the initial field distribution.Therefore a special to be iterated until input and output values of A coincide. The scalar formulation presented here avoids this restric- program module was included which calculates a rough (Gaussian-Hermite-Gaussian)approximation of the true tion. field distribution [24]from which the optimum choice of the artificial boundary and some design parameters of an C.Solution of the Algebraic Eigenvalue Problem optimized mesh division can be extracted.This program Equation (32)can be solved numerically using a suit- also generates a field distribution vector very well suited able standard algorithm.A good approximation of the true as an initial guess for the iterative procedure field function requires a high number of nodes.On the A typical grid as used for the finite-element calculations other hand,the overlap between two shape functions S is shown in Fig.4.For symmetry reasons the grid covers and S.is zero for nodes p and g belonging to different only a part of the half plane0.It is composed of four area elements.Therefore A and B are large,sparse ma- domains each of which consists of a number of rectangles. trices. Each rectangular element consists of two equally sized The coordinate-overrelaxation method [23]was used to triangular elements.The rectangle size within each of the four domains is compute the maximum effective indices and the associ- ated field distributions of the fundamental quasi-TE and domain A:4040 rectangles,0.2 x 0.1 gm quasi-TM modes of the waveguiding structures discussed domain B:505 rectangles,2.0 x 0.1 um above.This method completely exploits the sparseness of the matrices,thus saving memory and computation time domain C:400 rectangles,0.2 x 1.0 um By taking advantage of waveguide symmetry with respect to the y-axis and imposing a suitable boundary condition domain D:50 rectangles,2.0 x 1.0 um. for x =0,the quasi-TEo0.-TE1o.-TMoo.and -TMio Of course the grid can also be divided into a larger modes can be computed,where the subscripts designate number of domains with different element size.Thus the the mode order with respect to the width-and depth-co- mesh division can be optimized for the application to al- ordinates. most arbitrary waveguide geometries. In(20)the integration extends over the entire x-y plane Because guided mode fields of dielectric waveguides pos- V.NUMERICAL RESULTS sess an exponential decay for large distances from the By using the two procedures described in Sections III waveguide axis,we may introduce an artificial boundary and IV,we have computed and carefully compared the on which we impose the condition =0.If this boundary results obtained for the field distributions and for the ef- is chosen sufficiently far away from the waveguide axis, fective indices in a wide range of waveguide parameters. no significant error is introduced. As described in Section IV the results of the FEM calcu- The appropriate choice of the location of this boundary lations can be considered to be exact for the waveguides can be a serious problem because the maximum number examined in the present work.The comparison of the re- of area elements is limited by the computer memory size sults obtained by the ERI method with the FEM results A good geometrical resolution can be obtained using small therefore allows the verification of the approximations and elements;on the other hand,the total integration area also of some of the choices done with respect to the M. should be as large as possible in order to make the effect n.R.and k parameters of the almost analytical technique. of the artificial boundary as small as possible.This prob- Using the ERI technique we first approximated the re- lem was widely avoided by using small elements near the fractive index profile in depth and the effective index pro- waveguide axis,where the field distributions have large file in the transverse direction by (9)and (15)with M curvature.and large elements for the remaining region. R 1.With this choice the substrate indices of the actual

‘STRAKE ~f (11 GUIDED MODES OF TI LiNbO, CHANNEL WAVEGUIDES 1131 the propagation constant and field distributions of optical modes in the waveguide channel. An important fact is that A and B depend on the wave￾length X, but not on the effective refractive index Ne,,. Thus the effects of wavelength dispersion (cf., equation (6)) can be included very easily. This is an additional ad￾vantage of the scalar formulation as compared to vectorial techniques where, in general, the matrices depend on the effective refractive index while the wavelength is given by the eigenvalue. Then, in the vectorial formulation, due to wavelength dispersion of the refractive indices of the waveguide, the matrices will also depend on X, so Neff and X have to be specified as input parameters. As a result the eigenvalue will yield a wavelength value which, in gen￾eral, will be different from the input one. In order to ob￾tain a self-consistent vectorial solution the calculation has to be iterated until input and output values of h coincide. The scalar formulation presented here avoids this restric￾tion. C. Solution of the Algebraic Eigenvalue Problem Equation (32) can be solved numerically using a suit￾able standard algorithm. A good approximation of the true field function requires a high number of nodes. On the other hand, the overlap between two shape functions Sp and S, is zero for nodes p and q belonging to different area elements. Therefore A and B are large, sparse ma￾trices. The coordinate-overrelaxation method [23] was used to compute the maximum effective indices and the associ￾ated field distributions of the fundamental quasi-TE and quasi-TM modes of the waveguiding structures discussed above. This method completely exploits the sparseness of the matrices, thus saving memory and computation time. By taking advantage of waveguide symmetry with respect to the y-axis and imposing a suitable boundary condition for x = 0, the quasi-TEoo, -TElo, -TMoo, and -TMlo modes can be computed, where the subscripts designate the mode order with respect to the width- and depth-co￾ordinates. In (20) the integration extends over the entire x-y plane. Because guided mode fields of dielectric waveguides pos￾sess an exponential decay for large distances from the waveguide axis, we may introduce an artificial boundary on which we impose the condition 4 = 0. If this boundary is chosen sufficiently far away from the waveguide axis, no significant error is introduced. The appropriate choice of the location of this boundary can be a serious problem because the maximum number of area elements is limited by the computer memory size. A good geometrical resolution can be obtained using small elements; on the other hand, the total integration area should be as large as possible in order to make the effect of the artificial boundary as small as possible. This prob￾lem was widely avoided by using small elements near the waveguide axis, where the field distributions have large curvature, and large elements for the remaining region. Yt Fig. 4. Example of a typical grid used for the FEM computations. The coordinate-overrelaxation method is an iterative technique of solvfng (32); therefore, an initial guess for the field vector 4 is required. The number of iterations needed to achieve convergence strongly depends on the choice of the initial field distribution. Therefore a special program module was included which calculates a rough (Gaussian-Hermite-Gaussian) approximation of the true field distribution [24] from which the optimum choice of the artificial boundary and some design parameters of an optimized mesh division can be extracted. This program also generates a field distribution vector very well suited as an initial guess for the iterative procedure. A typical grid as used for the finite-element calculations is shown in Fig. 4. For symmetry reasons the grid covers only a part of the half plane x 2 0. It is composed of four domains each of which consists of a number of rectangles. Each rectangular element consists of two equally sized triangular elements. The rectangle size within each of the four domains is domain A: 4040 rectangles, 0.2 X 0.1 pm domain B: 505 rectangles, 2.0 X 0.1 pm domain C: 400 rectangles, 0.2 X 1.0 pm domain D: 50 rectangles, 2.0 x 1.0 pm. Of course the grid can also be divided into a larger number of domains with different element size. Thus the mesh division can be optimized for the application to al￾most arbitrary waveguide geometries. V. NUMERICAL RESULTS By using the two procedures described in Sections I11 and IV, we have computed and carefully compared the results obtained for the field distributions and for the ef￾fective indices in a wide range of waveguide parameters. As described in Section IV the results of the FEM calcu￾lations can be considered to be exact for the waveguides examined in the present work. The comparison of the re￾sults obtained by the ER1 method with the FEM results therefore allows the verification of the approximations and also of some of the choices done with respect to the M, 7, R, and K parameters of the almost analytical technique. Using the ER1 technique we first approximated the re￾fractive index profile in depth and the effective index pro￾file in the transverse direction by (9) and (15) with M = R = 1. With this choice the substrate indices of the actual

1132 JOURNAL OF LIGHTWAVE TECHNOLOGY.VOL,6.NO.6.JUNE 1988 TABLE III No.I cut polarization : (um] md V/un or A/un 19.66 0.417 4.8 2.9 6,536 (a) 0.6 4.135 3.50 4.714 18.17 0.385 4.9 2.9 6.532 ) 18.52 0.393 4.9 2.9 6.532 (c) 5.89 0.177 9.7 6.5 1.774E-2 (a) 多 1.2 4.135 4.135 3.331 5.42 0.162 11.5 6.3 1.70E-2 (b) 5.45 0.163 10.6 6.3 1.76E-2 (c) 1.90 0.055 9.7 6.8 1.691E-2 (a) 3 TH 1.2 4.135 4.135 3.429 1.93 0.056 12.9 6.6 1.51e-2 (b】 1.28 0.037 11.5 6.4 1.60E-2 (c) 0.58 0.014 9.2 7.0 2.846 (a) 1.2 4.135 3.50 4.051 1.35 0.033 14.1 6.6 2.53 (b) 0.37 0.009 11.8 6.2 2.87 (c) and the simulated waveguides coincide.The free param- approximation just explained (M =R 1);(c)will be eters n and k were chosen so as to achieve point matching discussed below.The reported examples range from well of the actual refractive index distribution with the approx- confined modes (b-values of about 0.4)to near cutoff (b imated one.We preferred to optimize the fitting in the nearly zero )It can be observed that the error in b.with central part (i.e.,near the z-axis,cf.Fig.1)of the re- respect to the FEM data,is negative for b>0.1 and fractive index distributions in order to minimize the dis- becomes positive for smaller b values.The behavior of crepancies for the well-confined modes.Point matching the error is closely related to the approximation of the in depth was imposed at u=1.0 and 1.5 for one and on. refractive index profile.In particular,if we require a good respectively.The "effective index profile''n(s)is sen- matching in the tail region of the on profile,an error arises sitive to the depth mode order and to the wavelength dis- in the central part of the curve (important for strongly persion of the field distribution in depth.Here,the match- confined low order modes);if we try to approximate well ing point was chosen at s =0.5sm,where s is taken from the central part,a large discrepancy occurs in the tail in- u'(s=Sm）+q=0 (33) troducing errors near cutoff. In Figs.5 to 8 the field distributions obtained by ERI which is the condition for mode cutoff with respect to the and FEM for the four cases considered in Table III,are depth coordinate. compared.The figures show sections through the field Some results obtained under these conditions are re- maximum.A very good coincidence is obtained for well ported in Table III together with the corresponding FEM guided modes;for modes near cutoff the main discrep- results.They refer to a waveguide obtained by diffusion ancies arise for the lateral field distribution X(x). of a Ti strip of width W=7 um and thickness 7 =274 In order to improve the overall accuracy of the ERI A before diffusion.The different cases of isotropic and technique we approximated the depth distribution of the anisotropic diffusion and modes of the two polarizations refractive index profile and the effective index profile at an operating temperature T 300K have been consid- nerr(s)in the range 0.1 6n(0.0)<on 6n(0.0)by ered. using values of M and R different from 1.This introduces Table III summarizes the following data:waveguide di- an error in the substrate index that invalidates the method mensions (D.,D )maximum refractive index change for modes near cutoff,but can allow a better representa- 6n(0,0)=n(0,0)-n,computed TEoo and TMoo mode tion for modes not too close to cutoff.The procedure we characteristics given as 6Nem=Nen-np,normalized used will now be presented. propagation constant b =(Ne-n)(n(0,0)-n) The choice of M and n used to approximate the function 6Nett/on(0,0).FWHM spot sizes of the field distri-exp (-Tu2)with r =I for one and I =Y for 6n.has butions in width (I)and in depth (I )maximum value been done definitively.After some trial the values M= of the normalized field strength max for P =1 W). 1.05 (for both cases),n 1.035 (T 1)andn 0.776 In Table III the reference solution (FEM)is labeled (a),(T =Y)were chosen.Such values correspond to point while the label (b)shows the results obtained with the matching the one and on,distributions with(9)at u=1.0

1132 JOURNAL OF LIGHTWAVE TECHNOLOGY. VOL. 6. NO. 6. JUNE 1988 and the simulated waveguides coincide. The free param￾eters r and K were chosen so as to achieve point matching of the actual refractive index distribution with the approx￾imated one. We preferred to optimize the fitting in the central part (i.e., near the z-axis, cf. Fig. 1) of the re￾fractive index distributions in order to minimize the dis￾crepancies for the well-confined modes. Point matching in depth was imposed at U = 1 .O and 1.5 for An, and 6n,, respectively. The “effective index profile” nerf( s) is sen￾sitive to the depth mode order and to the wavelength dis￾persion of the field distribution in depth. Here, the match￾ing point was chosen at s = 0.5s,,,, where s, is taken from p’(s = s,,,) + q’ = 0 (33) which is the condition for mode cutoff with respect to the depth coordinate. Some results obtained under these conditions are re￾ported in Table 111 together with the corresponding FEM results. They refer to a waveguide obtained by diffusion o! a Ti strip of width W = 7 pm and thickness 7 = 274 A before diffusion. The different cases of isotropic and anisotropic diffusion and modes of the two polarizations at an operating temperature T = 300K have been consid￾ered. Table I11 summarizes the following data: waveguide di￾mensions (Dr, D, ), maximum refractive index change 6n (0,O) = n (0,O) - nh, computed TE,, and TMoo mode characteristics given as &Ne, = Ne‘!, - n6, normalized propagation constant b = (Nzff - n,;)/(n2(0, 0) - n;) = 6Neff/6n (0, O), FWHM spot sizes of the field distri￾butions in width (r,) and in depth (r, ), maximum value of the normalized field strength +,,, ( for P = 1 W ). In Table 111 the reference solution (FEM) is labeled (a), while the label (b) shows the results obtained with the approximation just explained (M = R = 1 ); (c) will be discussed below. The reported examples range from well confined modes (b-values of about 0.4) to near cutoff (b nearly zero). It can be observed that the error in b, with respect to the FEM data, is negative for b > 0.1 and becomes positive for smaller b values. The behavior of the error is closely related to the approximation of the refractive index profile. In particular, if we require a good matching in the tail region of the 6n profile, an error arises in the central part of the curve (important for strongly confined low order modes); if we try to approximate well the central part, a large discrepancy occurs in the tail in￾troducing errors near cutoff. In Figs. 5 to 8 the field distributions obtained by ER1 and FEM for the four cases considered in Table 111, are compared. The figures show sections through the field maximum. A very good coincidence is obtained for well guided modes; for modes near cutoff the main discrep￾ancies arise for the lateral field distribution X(x). In order to improve the overall accuracy of the ER1 technique we approximated the depth distribution of the refractive index profile and the effective index profile neff(s) in the range 0.1 6n(0, 0) < Sn < 6rz(0, 0) by using values of M and R different from 1. This introduces an error in the substrate index that invalidates the method for modes near cutoff, but can allow a better representa￾tion for modes not too close to cutoff. The procedure we used will now be presented. The choice of M and q used to approximate the function exp ( -Tu2) with l7 = 1 for Sn, and r = y for An,,, has been done definitively. After some trial the values M = 1.05 (for both cases), 7 = 1.035 (r = 1 ) and 77 = 0.776 (I? = y) were chosen. Such values correspond to point matching the 6n, and 6n,, distributions with (9) at U = 1 .O

STRAKE er a/.:GUIDED MODES OF Ti:LiNbO,CHANNEL WAVEGUIDES 1133 100 100 80 60 60 0 20 40 20 20 12 u=15 00 04 08 245 81012 0 2 Fig.9.Typical behavior of n(s)(- -）and its approximation(···) -Y/um x /um used in the ERI method. (a) (b) Fig.5.Comparison between field distributions obtained using FEM (-) and ERI method with M=RI (and with M=1.05 (x)for waveguide 1 in Table 111.(a)Depth and (b)width field distribution. 07 100、 160 03 60 01 40 40 05 10 X/um Fig.10.Dispersion of the normalized effective refractive index for the 0 0 quasi-TMoo mode in the Y-cut waveguide having the following parame- 0 5 1015 20 25 30 0 5 10 tes:W=7um.D.4.1354m,D、=3.5m( -)andD、=4.I35 -Y/um x/um- um (---)FEM results are shown for D.3.5 gm (and D. (b) 4.135m(×). Fig.6.Comparison between field distributions for waveguide 2 in Table III. The values of R and k,used to approximate ne(s).de- pend on the wavelength and are chosen in the following 0■ 100 way:referring to the typical behavior of n(s)repre- sented in Fig.9,the program evaluates sm from(33);then 180 nento,R and k are computed by point matching (14)and 60, 60, (15)for s 0.s =5m/2 and s =(3/4)sm.With this * 0 is typically less than 2 percent.Usually M and R are u:15 larger than unity.Therefore we formally obtain"guided' oL 510152025 30 00 5 10 modes having NerI that produces a stronger guiding. 。。 20 Fig.10 shows b as a function of A for a TMoo mode in u-15 a Y-cut crystal.The agreement between the ERI-and the 0 20 30 0 0 10 20 FEM-results obviously increases as A decreases.The typ- x/um-- ical error in b ranges from nearly 50 percent for b=0.05 (a) (b) to approximately 4 percent for b =0.5. Fig.8.Comparison between field distributions for waveguide 4 in Table Also for the fields a significant discrepancy can be noted for b <0.1.The discrepancies in the field distributions occur in the evanescent tail;in depth,they arise always and 1.5,respectively.With these choices the discrepancy for u 1.5 (2.0)for the mode sensitive to on (on )i.e., between(9)and the Gaussian function is less than 3 per- in the regions where the approximation of on with the hy- cent for on 0.2 6n(0,0)and the error does not exceed perbolic tangent is not satisfactory.Even for modes near 20 percent forδn≥0.026n(0,0). cutoff the error in the depth spot-size is rather small (<8

STRAKE er U/.: GUIDED MODES OF Ti:LiNbO, CHANNEL WAVEGUIDES 20 1 I33 u=15 I I,, , 1 07 05 I D 03- 01 ; LO U“ {M 20 0 02L6 - xipm --c -Y/pm - (a) (b) Fig. 5. Comparison between field distributions obtained using FEM (-), and ER1 method with M = R = I (0) and with M = 1.05 (x) for waveguide 1 in Table Ill. (a) Depth and (b) width field distribution. 100 h 100 0 LO ! 20 U-2 0 5 10 15 20 25 30 -Ylum - xiurn - (a) (b) Fig. 6. Comparison between field distributions for waveguide 2 in Table I11 . 50’ x -13 loo 80 Fig. 7. Comparison between field distributions for waveguide 3 in Table 111. I 1 100 h loot 2 I 2x 20 Lo: \ \. I I I I 0 10 20 30 LO 50 0 10 20 30 -Y/um - x/pm - (a) (b) Fig. 8. Comparison between field distributions for waveguide 4 in Table 111. and 1.5, respectively. With these choices the discrepancy between (9) and the Gaussian function is less than 3 per￾cent for 6n > 0.2 6n (0, 0) and the error does not exceed 20 percent for 6n 2 0.02 6n(0, 0). 00 OL 08 12 Fig. 9. Typical behavior of &(s) (-) and its approximation (. . . ) used in the ER1 method. 5- The values of R and K, used to approximate n,?fr( s), de￾pend on the wavelength and are chosen in the following way: referring to the typical behavior of n&(s) repre￾sented in Fig. 9, the program evaluates s,, from (33); then n&,, R and K are computed by point matching (14) and (15) for s = 0, s = s,,/2 and s = (3/4) s,. With this choice the discrepancy in the range where n:rf( s) - ni > 0 is typically less than 2 percent. Usually M and R are larger than unity. Therefore we formally obtain “guided” modes having Ne,, 1 that produces a stronger guiding. Fig. 10 shows b as a function of X for a TMoo mode in a Y-cut crystal. The agreement between the ERI- and the FEM-results obviously increases as X decreases. The typ￾ical error in b ranges from nearly 50 percent forb = 0.05 to approximately 4 percent for b = 0.5. Also for the fields a significant discrepancy can be noted for b 1.5 (2.0)forthemodesensitiveto6n,(6n,,), i.e., in the regions where the approximation of 6n with the hy￾perbolic tangent is not satisfactory. Even for modes near cutoff the error in the depth spot-size is rather small ( < 8

1134 JOURNAL OF LIGHTWAVE TECHNOLOGY.VOL.6.NO.6.JUNE 1988 percent),and the method predicts the position and the 16]M.Fukuma.J.Noda,and H.Iwasaki.Optical properties of tita value of the field maximum very accurately.The discrep- nium-diffused LiNbO,strip waveguides,'J.Appl.Phys..vol.49. ancies arising in the lateral field distribution X(x)are pp.3693-3698.1978 17]G.P.Bava.I.Montrosset,W.Sohler,and H.Suche,'Numerical higher and usually a wider spot size is obtained.In fact modeling of Ti:LiNbO,integrated optical parametric oscillators." the solution of (8b)deals with a refractive index distri IEEE J.Quantum Electron..vol.QE-23,pp.42-51,1987. 8]M.V.Hobden and J.Warner.The temperature dependence of the bution n(s)becoming constant and equal to ng for s2 refractive indices of pure lithium niobate.Phys.Ler.vol.22.pp. sm while the actual index profile approaches n for s 243-244.1966. [9]H.Luidtke,W.Sohler,and H.Suche,"Characterization of Therefore the simulated waveguide used in the ERI Ti:LiNbO,optical waveguides,in Dig.Workshop Integrated Op- method has a smaller width than the original waveguide tics."R.Th.Kersten and R.Ulrich.Eds. Berlin,W.Germany as used by FEM and the lateral confinement is reduced. Technical Univ.Berlin,1980,pp.122-126. The choice of M and R values different from unity does 110]M.Minakata,S.Saito,and M.Shibata.'Two-dimensional distri bution of refractive index changes in Ti-diffused LiNbO,strip wave- not overcome this problem but allows better results for guides,J.Appl.Phys..vol.50.pp.3063-3067,1979 the normalized propagation constant b in the range of pa- [11]S.Fouchet,R.Guglielmi.A.Yi-Yan.and A.Carenco,Wavelength rameters where the ERI method is valid,i.e.,in case an dispersion of Ti:LiNbO:waveguides.'in Proc.2nd Eur.Conf.In- tegrated Optics (Florence.Italy).1983.pp.50-52. accurate approximation of the refractive index profile [12]G.Arvidsson and F.Laurell.Nonlinear wavelength conversion in mainly near the waveguide axis is required. Ti:LiNbO,waveguides,'Thin Solid Films.vol.136.p.29.1986. 113]T.Suhara.Y.Handa,H.Nishihara.J.Koyama,'Analysis of optical VI.CONCLUSION channel waveguides and directional couplers with graded-index pro file.''J.Opt.Soc.Amner.,vol.69,pp.807-815,1979. In this paper two methods for the calculation of optical [14]M.Koshiba.K.Hayata.and M.Suzuki.Approximatc scalar finite element analysis of anisotropic optical waveguides with off-diagonal modes of Ti:LiNbO:channel waveguides have been pre- elements in a permittivity tensor.IEEE Trans.Microwave Theory sented.The results computed using the purely numerical Tech.vol.MTT-32.pp.587-593.1984. finite-element method and the nearly analytical effective- 15]G.B.Hocker and W.K.Burns.Mode dispersion in diffused channel refractive-index method have been compared.The ap- waveguides by the effective index method,Appl.Opr.,vol.16.pp. 113-118.1977. proximations used in the ERI method lead to analytical 116]W.Streifer,D.R.Scifres,and R.D.Burnham.''Analysis of gain- expressions for the field distributions and effective in- induced waveguiding in stripe geometry diode lasers.'/EEE J. dices.The accuracy is generally good for well guided Quantum Electron..vol.QE-14.pp.418-427.1978. 117]H.A.Haus and R.V.Schmidt.''Approximate analysis of optical modes.The propagation constant and the position of the waveguide grating coupling coetficients.'Appl.Opt..vol.15,pp. fundamental mode field maximum are very accurate also 774-781.1976. for modes near cutoff,whereas the field distributions.es- [18]M.Abramovitz and I.Stegun,Handbook of Mathemarical Functions New York:Dover,1972,ch.22. pecially in the lateral (width)direction,become less ac- 119]D.Marcuse,Theory of Dielectric Optical Waveguides. New York curate for weakly guided modes.The proposed index pro- Academic 974 file approximation in the ERI method therefore leads to a 120]K.Hayata.H.Koshiba.M.Eguchi.and M.Suzuki."'Novel finite- clement formulation without any spurious solutions for dielectric very fast technique for those applications where a large waveguides.Electron.Lett..vol.22.pp.295-296.1986. number of computations have to be performed with the 121]E.Strake,"Berechnung und Messung von Feldverteilungen optischer Moden in dielektrischen Wellenleitern'',diploma thesis.Universitat- main interest lying on data like propagation constant, GH-Paderborn,Paderborn.W.Germany,1984. mode size or position of the mode field maximum.This [22]K.Hayata,M.Koshiba,and M.Suzuki,"Lateral mode analysis of case typically arises in the design and optimization of in- buried heterostructure diode lasers by the finite-element method.' tegrated optical devices. IEEE J.Quantum Electron..vol.QE-22.pp.781-788.1986. [23]H.R.Schwarz.Methode der finiten Elemente. Stuttgart.W.Ger- In all problems requiring an exact numerical solution of many:Teubner,1980. the wave equation while computation time is unimportant [24]S.K.Korotky,W.J.Minford.L.L.Buhl.M.D.Divino.and R.C. a numerical technique like the finite-element method must Alferness."Mode size and method for estimating the propagation constant of single-mode Ti:LiNbO:strip waveguides EEE be used.Its main advantage besides its high accuracy is Quantum Electron..vol.QE-18,pp.1796-1801,1982. that it can be applied to waveguides with arbitrary refrac- tive index profiles. Engelbert Strake was bor in Neuenkirchen/Westfalen.W.Germany,in REFERENCES 1958.He received the M.S.degree in physics from the University of Pa- derborn.Paderborn.W.Germany.in 1984. [1]R.A.Steinberg and T.G.Giallorenzi.''Modal fields of anisotropic Since 1986 he has been working for the Ph.D.degree at the University channel waveguides,"J.Opt.Soc.Amer..vol.67.pp.523-533. of Paderborn.Department of Applied Physics.in the field of integrated 1977. optics 12]C.Yeh,K.Ha,S.B.Dong,and W.P.Brown.'Single-mode optical waveguides.App/.Opr..vol.18.pp.1490-1504.1979. 13]N.Mabaya,P.E.Lagasse.and P.Vandenbulcke,Finite element analysis of optical waveguides,'IEEE Trans.Microwave Theory Gian Paolo Bava was born in Varallo,Italy,in 1937.He received the Trh.vol.MTT-29,pp.600-605.1981. degrec in electrical engineering from Politecnico di Torino.Torino.Italy. 14]J.-D.Decotignie,O.Parriaux.and F.E.Gardiol,Wave propaga in1961, tion in an inhomogeneous diffused channel guide using a finite-differ. Since 1961 he has been with the Institute of Electrical Communications ence technique,Proc.Ist Eur.Conf.Integrated Opt.(London,En- (succesively Department of Electronics),Politecnico di Torino.From 1970 gland).1981.Pp.40-42, to 1976 he was in charge of the course on microwave techniques at the 15]W.K.Burns,D.H.Klcin.and E.J.West.Ti diffusion in Department of Electronics.At present his main research activities are con- Ti:LiNbO,planar and channel optical waveguides.''J.Appl.Phys.. cerned with microwave MESFET amplifiers and oscillators and integrated vol.50,pp.6175-6182.1979 optics

1134 JOURNAL OF LIGHTWAVE TECHNOLOGY. VOL. 6. NO 6, JUNE 1988 percent), and the method predicts the position and the value of the field maximum very accurately. The discrep￾ancies arising in the lateral field distribution X(x) are higher and usually a wider spot size is obtained. In fact the solution of (8b) deals with a refractive index distri￾bution n&( s) becoming constant and equal to ni for s I s,,~ while the actual index profile approaches nb for s -+ W. Therefore the simulated waveguide used in the ER1 method has a smaller width than the original waveguide as used by FEM and the lateral confinement is reduced. The choice of M and R values different from unity does not overcome this problem but allows better results for the normalized propagation constant b in the range of pa￾rameters where the ER1 method is valid, i.e., in case an accurate approximation of the refractive index profile mainly near the waveguide axis is required. VI. CONCLUSION In this paper two methods for the calculation of optical modes of Ti : LiNbO, channel waveguides have been pre￾sented. The results computed using the purely numerical finite-element method and the nearly analytical effective￾refractive-index method have been compared. The ap￾proximations used in the ER1 method lead to analytical expressions for the field distributions and effective in￾dices. The accuracy is generally good for well guided modes. The propagation constant and the position of the fundamental mode field maximum are very accurate also for modes near cutoff, whereas the field distributions, es￾pecially in the lateral (width) direction, become less ac￾curate for weakly guided modes. The proposed index pro￾file approximation in the ER1 method therefore leads to a very fast technique for those applications where a large number of computations have to be performed with the main interest lying on data like propagation constant, mode size or position of the mode field maximum. This case typically arises in the design and optimization of in￾tegrated optical devices. In all problems requiring an exact numerical solution of the wave equation while computation time is unimportant a numerical technique like the finite-element method must be used. Its main advantage besides its high accuracy is that it can be applied to waveguides with arbitrary refrac￾tive index profiles. [I1 121 131 141 151 REFERENCES R. A. Steinberg and T. G. Giallorenzi. “Modal fields of anisotropic channel waveguides,” J. Opt. Soc. Amer., vol. 67, pp. 523-533, 1977. C. Yeh, K. Ha, S. B. Dong, and W. P. Brown, “Single-mode optical waveguides,” Appl. Opt., vol. 18, pp. 1490-1504, 1979. N. Mabaya, P. E. Lagasse, and P. Vandenbulcke, “Finite element analysis of optical waveguides,” fEEE Trans. Microwave Theory Tech., vol. MTT-29, pp. 600-605, 1981. J:D. Decotignie, 0. Parriaux, and F. E. Gardiol, “Wave propaga￾tion in an inhomogeneous diffused channel guide using a finite-differ￾ence technique,” Proc. 1st Eur. Con$ fntegrated Opt. (London, En￾gland), 1981, pp. 40-42. W. K. Burns, D. H. Klein, and E. J. West, “Ti diffusion in Ti : LiNbO, planar and channel optical waveguides,” J. Appl. Phys., vol. 50, pp. 617556182. 1979. M. Fukuma, J. Noda. and H. lwasaki, “Optical properties of tita￾nium-diffused LiNbO, strip waveguides.” 1. Appl. Phys.. vol. 49, pp. 3693-3698, 1978. G. P. Bava, I. Montrosset. W. Sohler. and H. Suche, “Numerical modeling of Ti : LiNbO, integrated optical parametric oscillators,” fEEEJ. Quantum Electron.. vol. QE-23, pp. 42-51, 1987. M. V. Hobden and J. Warner, “The temperature dependence of the refractive indices of pure lithium niobate,” Phvs. Lett., vol. 22, pp. H. Ludtke, W. Sohler, and H. Suche, “Characterization of Ti : LiNbO? optical waveguides,” in Dig. Workshop Integrared Op￾tics,” R. Th. Kersten and R. Ulrich, Eds. Berlin, W. Germany: Technical Univ. Berlin, 1980, pp. 122-126. M. Minakata, S. Saito, and M. Shibata, “Two-dimensional distri￾bution of refractive index changes in Ti-diffused LiNbO, strip wave￾guides,” .I. Appl. Phys., vol. 50, pp. 3063-3067, 1979. S. Fouchet, R. Guglielmi, A. Yi-Yan. and A. Carenco, “Wavelength dispersion of Ti:LiNbO, waveguides,” in Proc. 2nd Eur. Con$ In￾tegrated Optics (Florence, Italy), 1983, pp. 50-52. G. Arvidsson and F. Laurell. “Nonlinear wavelength conversion in Ti: LiNbO, waveguides,” Thin Solid Films, vol. 136. p. 29, 1986. T. Suhara, Y. Handa, H. Nishihara, J. Koyama, “Analysis of optical channel waveguides and directional couplers with graded-index pro￾file.” J. Opt. Soc. Amer., vol. 69, pp, 8077815, 1979. M. Koshiba, K. Hayata, and M. Suzuki, “Approximate scalar finite￾element analysis of anisotropic optical waveguides with off-diagonal elements in a permittivity tensor,“ fEEE Truris. Microwave Theory Tech.. vol. MTT-32. pp. 587-593, 1984. G. B. Hocker and W. K. Burns, “Mode dispersion in diffused channel waveguides by the effective index method,” Appl. Opt., vol. 16, pp. W. Streifer, D. R. Scifres. and R. D. Burnham. “Analysis of gain￾induced waveguiding in stripe geometry diode lasers.” fEEE J. Quuntum ELectron., vol. QE-14, pp. 418-427, 1978. H. A. Haus and R. V. Schmidt, “Approximate analysis of optical waveguide grating coupling coefficients,” Appl. Opt., vol. 15, pp. M. Abrdmovitz and I. Stegun, Hundhook ofMuthc~n~aticul Functions. New York: Dover, 1972, ch. 22. D. Marcuse, Theory of Dielectric OpticNI Wuvrguidrs. New York: Academic, 1974. K. Hayata, H. Koshiba. M. Eguchi. and M. Suzuki, “Novel finite￾element formulation without any spurious solutions for dielectric waveguides,” Electron. Letr., vol. 22, pp. 295-296, 1986. E. Strake, “Berechnung und Messung von Feldverteilungen optischer Moden in dielektrischen Wellenleitern”. diploma thesis, Universitit￾GH-Paderborn, Paderborn, W. Germany, 1984. K. Hayata. M. Koshiba, and M. Suzuki, “Lateral mode analysis of buried heterostructure diode lasers by the finite-element method.” IEEE J. Quanrum Electrori.. vol. QE-22, pp. 781-788, 1986. H. R. Schwarz. Merhork drrjitiircw Elemente. Stuttgart. W. Ger￾many: Teubner, 1980. S. K. Korotky, W. J. Minford. L. L. Buhl, M. D. Divino, and R. C. Alferness. “Mode size and method for estimating the propagation 243-244, 1966. 113-118, 1977. 774-781, 1976. constant of single-mode TI : LiNbO? strip waveguides,.” iEEE J. Qucnitum Electron., vol. QE- 18, pp. 1796- I 80 1, 1982. * Engelbert Strake was born in NeuenkircheniWestfalen, W. Germany, in 1958. He received the M.S. degree in physics from the University of Pa￾derborn, Paderborn, W. Germany, in 1984. Since 1986 he has been working for the Ph.D. degree at the University of Paderborn, Department of Applied Physics, in the field of integrated optics. * Gian Paulo Bava was born in Varallo, Italy, in 1937. He received the degree in electrical engineering from Politecnico di Torino, Torino, Italy, in 1961. Since 1961 he has been with the lnstitute of Electrical Communications (succesively Department of Electronics), Politecnico di Torino. From 1970 to 1976 he was in charge of the course on microwave techniques at the Department of Electronics. At present his main research activities are con￾cerned with microwave MESFET amplifiers and oscillators and integrated optics

STRAKE a/.:GUIDED MODES OF Ti:LiNbO,CHANNEL WAVEGUIDES 1135 Ivo Montrosset was bom in Aosta.Italy,in 1946.He received the degree 1972 to 1986.In 1982 he became an Associate Professor.In June 1986 he in electronic engineering from the Politecnico di Torino,Torino,Italy,in was appointed Full Professor at the Department of Biophysics and Elec- 1071 tronics.University of Geneva,Geneva,Italy.His main activities are in the He was with the Department of Electronics.Politecnico di Torino.from field of numerical methods for wave propagation and guided-wave optics

STRAKE PI ul.: GUIDED MODES OF Ti: LiNbO, CHANNEL WAVEGUIDES 1 I35 Ivn Montrosset was born in Aosta, Italy, in 1946. He received the degree in electronic engineering from the Politecnico di Torino, Torino, Italy, in 1971. He was with the Department of Electronics. Politecnico di Torino, from 1972 to 1986. In 1982 he became an Associate Professor. In June 1986 he was appointed Full Professor at the Department of Biophysics and Elec￾tronics, University of Geneva, Geneva, Italy. His main activities arc in the field of numerical methods for wave propagation and guided-wave optics