- Email: [email protected]
Analysis of optical fibers with arbitrary refractive index profiles: accuracy, convergence, and effects of finite cladding
Optics Communications North-Holland
99 ( 1993) 393-404
OPTICS
COMMUNICATIONS
Full length article
Analysis of optical fibers with arbitrary refractive index profiles: accuracy, convergence, and effects of finite cladding Lakshman
S. Tamil
Erik Jonsson School ofEngineering Richardson, TX 75080. USA
And Computer Science and
Centerfor Applied Optics, The University of Texas at Dallas,
and Gregory H. AickIen Erik Jonsson School of Engineering And Computer Science, The University of Texas at Dallas, Richardson, TX 75080, USA Received
24 September
1992; revised manuscript
received 4 January
1993
We have formulated a matrix eigenvalue problem for cylindrical optical fibers from a set of finite difference equations. Numerical solution of this problem yields the propagation constants for propagating modes. The method can be used for arbitrary index profiles, does not require the explicit evaluation of Bessel or modified Bessel functions, and does not use iterative methods to search for the propagation constants as was the case in earlier proposed methods using finite differences. The method is accurate, fast, and simple. We have established the convergence and stability of this method, and explored the effects of finite cladding width on the dispersion characteristics.
1. Introduction Wave propagation in optical fibers has been analyzed using various methods. We will be using a finite difference method. Other methods proposed to find the propagation constants of guided modes in optical fibers with arbitrary refractive index profiles include the WKBJ method, variational method, power series expansion method, staircase approximation method, and finite element method. The WKBJ method [ 1,2] is a geometrical optics approximation that works whenever the refractive index of fiber varies only slightly over distances of the order of the optical wavelength and are applicable only to thick fibers in which many modes can propagate. For those fibers in which only a few modes propagate, the error of the WKBJ method increases intolerably and this method is not applicable to modes near cutoff. Besides, the effect of an index valley at the core-cladding boundary, which plays an important role in reducing multimode dispersion, cannot be treated by the ordinary WKBJ method. In 0030-4018/93/%06.00
0 1993 Elsevier Science Publishers
the variational method the scalar wave equation is converted into a variational problem subject to the given boundary conditions. The variational problem is solved either by using the Rayleigh-Ritz method [ 31 or perturbation method [4]. In the RayleighRitz method the eigenfunction is expressed in terms of a set of orthogonal functions and the variational function is minimized. The disadvantage is that we need to assume a trial function [ 5 1. In the perturbation method of analysis, the computation of the propagation characteristics for an arbitrary profile is done by correcting the solution for a uniform core fiber considering the difference in the profile as the perturbation term. The power series expansion method [ 61 consists of expressing the refractive index for the field term by term. This method is useful only for cases in which the refractive index profile can be expressed by a relatively simple power series. In some cases the series do not converge and this method is not applicable [7]. In the staircase approximation method [8,9] the refractive index is approximated by an appro-
B.V. All rights reserved.
393
OPTICS
COMMUNICATIONS
staircase function. The ave equation is solved In each stratified layer and te solutions are then ‘.:onnccted at the cylindrical bt ndaries between these layers to obtain the proper sc ltion representing the propagation characteristics. he number of layers should be infinite in order tl .t the refractive index profile approaches that of act al fiber profiles. Thus, the results of the propagation instant will differ from the actual values when a fin e number of layers is used. A large number of laye requires considerable computer time and hence in .his method the accuracy and computer time are -aded off. The fiber problem has bee analyzed by Okamato and Okoshi [ lo] using a fin : element method formulated in the axial fields. ‘he problem with this method is that it suffers fror spurious modes when the finite elements are not arefully chosen [ 111. Lenahan [ 121 has formulat 1 a matrix eigenvalue problem from a finite elen :nt analysis using the Galerkin weighted residua method. To achieve computational efficiency, a 1 ecewise linear approximation to the solution func on must be used. In this paper we present 2 efficient finite difference method to find the proI gation constants of optical fibers with arbitrary rt ractive index profiles. The method does not involv a search procedure to find the propagation consta ts, nor does it require explicitly evaluating Besse and modified Bessel functions, as was the case in .he earlier works on finite difference analysis of 01 icalfibers [ 13,141. We construct a matrix equatior from a set of simultaneous finite difference ec .ations governing the propagation in an optical fi :r and solve for the eigenvalues to obtain the pr pagation constants. In sect. 2 we give the mathem; ical formulation of the discretized differential equat in at various grid points in the radial direction and t : construction of a matrix equation incorporating le boundary conditions at the core-cladding inteif :e and the jacket. Extending our method, which s formulated for a-index fibers, to arbitrary ref ctive index profiles is covered in sect. 3. In sect. 4 e discuss the numerical evaluation of propagation c nstants and present results. This includes a discu ion of the convergence and stability of the method ilong with the effect of the number of grid points I Lthe computation, and the effects of finite claddi g width on dispersion 1-n-k tr’
394
characteristics.
2. Mathematical
15 June 1993
Our conclusions
are given in sect. 5.
formulation
The optical fibers considered here are inhomogeneous dielectric cylinders of radius Q called the “core” surrounded by a homogeneous refractive index medium called the “cladding”. The cladding, in turn, is encased in a highly lossy material called the jacket. A representative fiber cross-section is shown in fig. 1. The refractive index profile of the fiber, called an a-index profile, is given by n2(1)=n:[1-2pd(r/a)“], =n:[1-2d])
forOa .
(1)
Here, n, is the maximum refractive index of the core, d is the relative refractive index difference between the core axis and cladding, p a parameter representing the refractive index step or valley at the corecladding boundary. A smooth continuation at the core-cladding boundary, the presence of a step, and that of a valley are expressed by p = 1, p < 1, and p > 1, respectively. {a 1CxeR} is a profile parameter. Some examples of I-u-index profiles are shown in fig. 2. The propagation characteristics of an optical fiber are governed by the scalar Helmholtz differential equation [ 15 ]
Jacket -L Fig. 1. Optical fiber showing grid points used in the example.
Volume 99, number
5,6
OPTICS
Core
15 June 1993
COMMUNICATIONS
Cladding
Region
1 .52
Step
rho Parabolic
Region
Index
1
=
Index
1.o Radius, r Fig. 2. Examples of a-index profiles. LY=co yields a step index, while (Y= 2 yields the parabolic index. Values of p control the characteristics of the interface between the core and the cladding; p< 1 results in a step at the interface, p> 1 yields a valley.
n2(r)k2-p2-
$
y=o.
This scalar wave equation is the simplification of the exact vector wave equation under the assumption that Vn/n is small, which includes the “small index gradient” and “weakly guiding approximations” [ 16,171. In the above equation v(r) is the transverse field function which may denote either the dielectric field or the magnetic field, r is the radial coordinate, n(r) is the radial refractive index profile, k is the vacuum wave number, fi is the propagation constant which is to be computed, and m is a mode parameter given by m=l,
forTEandTMmodes
(n=O),
=n+l,
for EH modes (n~n\l) ,
=n-1,
for HE modes (noN)
.
=O,
(2) Y(O)=0
We need to solve the differential equation in order to compute the propagation constants. From the rotational properties of w the associated boundary condition at the center of the core (r=O) is
formp0.
9
(4)
The other boundary condition applied tinction of field at the jacket written as vjacket
=
vru,= b =
0
is the ex-
>
(5)
where b is the radius of core and cladding
together.
2.1. Transformation to nondimensional form We need to nondimensionalize the differential equation for easy computation. This is achieved by setting u=wlv0,
(3)
form=O,
x=rla,
(6)
where v/Ois the maximum field amplitude and a is the radius of the core. Substituting eq. (6) into eq. (2) we obtain kzn2(xa)-/12-$$
>
u=o.
(7)
395
15 June 1993
OPTICS COMMUNICATIONS
by including the refractive i dex distribution given ly cq. ( I), the above equatil 1 can be rewritten as
( 11 ), and defining
-!Fy
CI’U
+a2 k2nT[ l-_/j
---; + C-Y-
h is the width between the grid points and x=ih, {i=O, 1, 2, . ..I. Substitutingeqs. (13) and (14) into
:a)]-/32-
?J
g
(
>
V2nT
17,
=: 0,
(15)
n:--nz’
(8)
we get
where
J‘(r)=2pA(r/a)U, =2A,
a
Defining
the parameters
,li=a(k2n:
-P2)
we can define
u~+l-2ut+ui_I + fUi+l-Ul_I
OQr
“2 )
h2
ed frequency,
propagatic
(9) I constant,
U-i=0 3
and on rearranging, as [ 171
V2= U2+ W2=k2a2(nT -ni and the modified
2h
+
C and Was
WE& P2-k2n:)‘/2,
V, the normali
rh
8, as
the equation
(16)
becomes
u,_i[-$(l-t>] +u,[$(2+ $) +Bliiha)-ij]
=a2(k2n:-/3’). +.,+,[-;(1+;)]=0.
(10) ‘Then eq. (8) becomes
p---V2n:
!e+;g+
(
?Z: -n:J
2 ax,->
u=o, 1
(11)
with
O,
f(x)=2pdxa,
(12)
x>l.
=2A,
2.2. Discretizing the differer; ial equation When the function u and is derivative are single functions of x, the first valued, finite and continuou an be approximated by and the second differentials third order difference form. las as follows [ 18 ] :
du
-z
2h
dx
d2u s=
(13)
’
h2
396
ui+,
+I?f(lha)-B)
+u&$
=O,
(18) where
=O,
m=O, m#O.
(14)
.
u* )
4+2m2-6 2h2
Ul
At i=2,
Ui+i -22ui+U,-i
Here ui=u(x)
For the purpose of illustration, we have chosen six grid points along the radial direction as shown in fig. 1. In general, the number of grid points can be any number not less than four, the minimum necessary to take care of the boundary conditions. Depending on whether m = 0 or m # 0, the field or its derivative vanishes at the center of the core. When the derivative of the field vanishes, u0 = u,. Writing finite difference equations at the grid points, we obtain the following set of equations. At i= 1,
6=1,
%+I -u,-I
(17)
=u(x+h)
u,_, =u(x-h)
,
(3 +u,($ =o.
+Ff(2ha)-j?)
+u&$ (19)
Volume 99, number
$6
OPTICS
COMMUNICATIONS
At i=3,
15 June 1993
- (2i+ 1) ai,,+
+Yf(3ha4
+u,($)
=o.
(20)
1 =
2ih2
(25)
’
In order to convert the problem into an eigenvalue problem, we rewrite eq. (24) as (T-fl)u=O,
At i=4,
=o.
(21)
At i=5, since the field goes to zero at the jacket,
(22) Finally, at i= 6, again using the boundary that the field goes to zero at the jacket,
(26)
where I is the identity matrix, and T is a tri-diagonal matrix. Equation (26 ) defines an eigenvalue problem. This means that eq. (26) has a nontrivial solution if and only if p”is an eigenvalue [ 191. Hence, the required normalized propagation constants contained in Fare obtained by finding the eigenvalues of the tri-diagonal matrix T. This mathematical formulation can be generalized to {n 1n 1EN} grid points in the radial direction of the fiber without difficulty.
condition
3. Arbitrary profiles, multipie layers, and field ug=o.
(23)
Since the boundary condition in eq. (23) is incorporated in eq. (22), we have a system of five equations. 2.3. Matrix equation formulation Formulating a matrix equation from the above set of equations for the convenience of generalization and easy computation, we obtain a,, -s a21
a12 a22
Au=
Ul
-P
a32
a23 a33
u2
-P
a43
a34 a44
-
u3 B
a54
i
a45 a55-B
=o,
iii
~5
(24)
where the matrix elements
ai+, = a,,i =
u4
- (2i-
are defined
by
1)
2ih2
2i2-!-2m2-6 2i2h2
’
+pf(iha),
2i2+m2 = i~h2 + pf(iha)
,
i=l
i#l,
,
distributions We have developed this method of analyzing optical fibers using the a-index profile. This is because the a-index profile is commonly used in the literature and can represent a large variety of real refractive index profiles, including the very important step and parabolic profiles. But our formulation is not limited to a-index profiles. To see how to extend the method to arbitrary profiles without rederiving a system of finite difference equations, consider eq. ( 11). The refractive index profile is included in this equation through the functionf(x), which is defined in eq. (12). Usingf(x), the refractive index profile of the fiber can be rewritten as
n2(x)=n:[
1 -f(x)]
Solving for f(x) f(x)=
1 -n2(x)/nf.
.
(27)
yields (28)
By generating the discretized A=f(x,= ih ) from samples of an arbitrary refractive index profile n(x,), the method we have outlined in this paper can be used directly on arbitrary profiles, as long as the “weakly guiding” approximation holds. Multiple layer waveguides of any number of layers may be considered special cases of arbitrary refrac-
397
‘~OlUllli 0
1’.
nurlbcr 5.‘)
OPTICS COMMUNICATIONS
t ive inde!x profiles Since WC’II Lve normalized the fiber c’#:)rcradius. ,z, to unity, ‘rv~’must explicitly define which layer:, comprise the )re before using our methcd and scale all quantit s accordingly. Since the propagation const mt for a mode i, /I$, is uniquely associated with a fi .he field distributions for a propag.sting mode can k determined from the cigenve’ctors u of matrix T ( I( e eq. (26) ). Many eigenvalue routines will return (’ genvectors as well, but at the cost of greatly incr’: lsing the number of computations. Useful approximations tcl the eigenvectors for propagating modes can be co. lputed by constructing the tri-diagonal matrix T, SI btracting a specific fl from each element of the ma:’ Ldiagonal, and solving for the melem:nts of u using st;l ldard techniques from linear algebra. From the ob:l rvation that for propagating modes the field will approach zero at the cladding/jacket boundary, V’ can set uN, the rightmost element of u, to a very ! mall value (not zero), and use a simple backsubst ution process to solve for the rest of the u,. This 11ocedure yields a good approximatj on to the field di! tribution multiplied by an arbi trar), constant.
4. Nwneriad
evaluation,
rest Its, and discussion
Although the mathemati: al formulation of our method for determining the propagation characteristics ad an optical fiber is co Lched in terms of matrix equations, there are special structures that lead to very efficient numerical imp mentations. First, since T is a tri-diagonal matrix, vi : can use sparse matrix techni’ques to reduce storal e requirements for T. Second, since T is a quas symmetric tri-diagonal matrix, we can use a sirnil lrity transformation to convert T into a real, symt metric matrix [20]. Finally, ,the ei genvalues of a re: I , symmetric matrix may be computed using an efficil: It 0 (N2 ) algorithm (in our case, the tq1i.c routine f [ ,rn ref. [ 2 11, which has an alperation count of apprr ximately 30N’). Using eqs. (25), we havt implemented a pair of C language: programs which :ompute the normalized propagatic’n constants for ! ibers with arbitrary refractive index profiles. WI: define the normalized propagation constant as 398
x_
;z _
15 June 1993
k2n?-D2 k*n: -k*n:
’
Note that some authors normalized propagation b=1-U2/V2=1-x.
(29) (e.g. Gloge [ 221) define a constant as (30)
One program computes x for all propagating modes at a specific value of m in eq. (2) over a range of normalized frequencies V. Another program searches for the cutoff frequency ( V,) of a specific linearly polarized (LP) mode. Both programs allow all of the parameters in eq. ( 1) to be varied, as well as the values of b and N,, which are the fiber radius (see fig. 1) and number of grid points in the fiber core, respectively. In verifying the performance of our method, we have computed the propagation characteristics of step index and parabolic index fibers over a normalized frequency range of 0 to 20. These index profiles have analytical solutions and have been studied analytically and numerically by other authors [ 14,23-261. Our results agree well with previously published results, as shown in table 1. Note that for propagating modes, x must lie between 0 and 1 (i.e. O
Volume 99, number
56
OPTICS
COMMUNICATIONS
15 June 1993
Table I Comparison of the cutoff frequencies obtained by the finite difference method with analytical and previous numerical results. Raber is the fiber radius, R,, is the core radius, and 6 is the percentage difference from the infinite cladding result. 256 points were used in the fiber core. cr
1 2
Infinite cladding
Mode (m, I)
4.381 3.518 7.45 1 5.744 9.645 7.848 9.904 3.181 3.000 2.886 2.649 2.527 2.405
1,l I,1 12
2,1 22
3 4 5 10 20 co
Normalized
3,l 4,l 1,l 1,l 171 1,l 1,l 121
1
cutoff frequency
R m-x= lOR,re
6 (%)
Rfiber
4.391 3.526 7.457 5.744 9.645 7.848 9.904 3.189 3.007 2.894 2.657 2.535 2.413
0.23 0.23 0.08 c lo-* < 1o-2 <10-z < 1o-2 0.3 0.2 0.28 0.30 0.32 0.33
4.384 3.520 7.453 5.744 9.645 7.848 9.904 3.183 3.002 2.888 2.65 1 2.529 2.407
=
2’%0,
6 (%) 0.07 0.06 0.03
.o
5
Fig. 3. Dispersion
characteristics
15
10
Normalized
Frequency,
of a step index fiber (a=a,
frequency using from 4 to 256 points in the core, and for fiber radii from 5 to 20 times the core radius. From this figure we can see the expected convergence on a final result as the number of points in the
20
V
A=O.O38). R,+._=IOR,,.
core increases. The effect of cladding width is also apparent. The effect of the number of grid points in the core is two fold. As the number of grid points is increased the distance between samples of the refractive index
399
‘L’ol~.lmc9’7, number 5,b
OPTICS
COMMUNICATIONS
5
10
Normalized Fig. 4. Dispersior
0.0
characteristics
c
’
’
’
I
I
I
I
I
5
400
characteristics
V
I,
I
I
10
Normalized Fig. 5. Dispersion
of a parabolic
20
15
Frequency,
index fiber (cu=2, p= 1,4=0.038).
ofa parabolic
’
15 June 1993
Frequency,
I
I 15
I
I
Rfikr= IOR,,.
!
I
I 20
V
index fiber (cy=2, p=O.75, d=O.O38). RBkr= IOR,,.
Volume 99, number
5,6
OPTICS
COMMUNICATIONS
15 June 1993
L
0.0
0
5
10
Normalized Fig. 6. Dispersion
characteristics
’ ’
-0.05 0
’
’
’ 50
ofa parabolic
’
’ ’ ’
’
100
’
15
Frequency,
20
V
index fiber ((r=2,p=2,4=0.038).
’ ’ ’
’ 150
’
’ ’ ’
’ 200
’ ’ ’
Rfibcr= lOR,.
’
’ 250
’
’
’ ’ 300
No. of Grid Potnts in the Core Fig. 7. Convergence behavior of computed cutoff frequency with the number of grid points in the fiber core for the LP,, and LPO* modes of a step index fiber. 401
Volume 99. number 5,6
OPTICS COMMUNICATIONS
profile is decreased, resulting in a better approximation of the actual profile. This is especially apparent in profiles with sharp t ansitions at the corecladding interface, such as fc r a!>> 1. Cases where p # 1 (see eq. ( 1) ) are also lil ely to be poorly modeled by a small number of co *e sample points. The effect of reducing the number if sample points in the core is to apply a “low pass” spatial filter to the refractive index profile. Setting the number of poir ts in the core also effectively applies a filter to the spatial frequency content of the field distributior s calculated for each mode in the fiber. When ccImputing propagation constants at higher norrnalizc d frequencies, using a small number of samples ma t induce errors due to a form of “aliasing”. These 1wo effects are responsible for the poor results wh :n the number of grid points in the core is below ap proximately 10 for the modes we have examined. Using the step index fiber as an example, fig. 8 demonstrates the effects of th : number of grid points by plotting the computed cute ff frequency for modes LPO,, LP, ,, and LP,, for sew ral different grid sizes. In this figure, each curve is r ormalized to the value
100
15 June 1993
of the cutoff frequency for that mode calculated with 256 points in the core. We can see that for mode LPol the cutoff frequency calculated with 8 points in the core is less than 1.0015 times that computed using 256 points in the core, while for mode LPI, (with a higher cutoff frequency) we need at least 12 points in the core for similar results. In general, as the normalized frequency increases, the number of points in the core must be increased to maintain the accuracy of the method. For modes with relatively low cutoff frequencies, variations in cladding width produce large changes in the calculated cutoff frequency, V,. Cutoff frequency increases as the cladding width decreases. This is the expected behavior. Analyses assuming infinite cladding width, while adequate for many purposes, fail to account for the increasing importance of finite cladding width as the normalized frequency becomes smaller. The fundamental mode, which has no cutoff frequency when the cladding is infinite, shows a definite cutoff in real fiber. Figure 9 shows the effects of cladding width on the cutoff frequencies of two LP modes in a step index fiber. In this plot, the curves for each mode are nor-
150
200
No. of Grid Points in the Core
Fig. 8. Effect of the number of gri i points on the computed cutoff frequencies of propagating modes of a step index fiber (LPO,, LPI,, and LP,,) .
402
15 June 1993
OPTICS COMMUNICATIONS
Volume 99, number 5,6
Cladding Width/Core
Width
Fig. 9. Effect of cladding width on cutoff frequency for a step index fiber.
0.0
0.5
1.5
1 .o
Radius,
2.0
I
Fig. 10. Field patterns for the LP 01, LP,I, and LPtl modes of a step index fiber with 4~0.038. These patterns were computed using R tibcr=1ORc.m V= 10,and 32 points in the core.
‘I!) lr:, : j’i. number
5.6
OPTICS
COMMUNICATIONS
#‘II:1):’c 8.110 the cutoff frequent for that mode in the nl 1.1‘1t .: cladding case. As an ex lmple, consider mode 1.J , , 14/t:can see that the cum ‘ffrequency when the I‘ilzer T.rdius is 10 times the COI : radius is only about I .(rO ; imes (0.3%) greater thz t the cutoff frequency wt. er I.hazfiber radius is 25 t mes the core radius. How e #Tel-,when the fiber radi ts is only 5 times the core mdius then the cutoff f equency increases to 1.O14 I.imes ( 1.4%) greater th: n the cutoff frequency when Ime fiber radius is 25 ti les the core radius. Usi’rg the approach outline d in sect. 3, we have camp .ned the field distributic IS for three LP modes in a !il repindex fiber. Figure 10 shows the results, with the -fis:ld patterns normalized so that the maximum VElue in each pattern is one. 1 hese results agree well u’ th, I he results reported in r-1f. [ 15 1.
5. Co ,dlusions We know that when d is sn all in an optical fiber, tit e s( al;ar approximation yiel Is results that are very close ‘to the exact vector forrr ulation. Even for large dlIfft:r lences between the core .nd cladding refractive itIde> es, optimum single-mot e fiber parameters obt;.inel:I ffrom the scalar apprc ximation differ negligjbly from those obtained u! _ng the exact formulation [ 271. V& heave developed a methl d to evaluate the propagatii:m constants by transfc .ming the scalar wave equa ;ion into a set of finite d fference equations and then converting into a matr x eigenvalue problem. The imethod does not involv a search procedure to l‘nd I.he propagation constat ts, or the explicit evalIIatic n of Bessel and modified Bessel functions, which i:s tirae consuming, as was tf : case in earlier works. W;: lhave demonstrated t te convergence of the tnetl~od and the dependence of the rate of conver,g;ence on the number of grit points. The method is eleg; nt., stable, straight forw .trd, is applicable to arbitrary index profiles, and ir accurate. We have also expllored and established the effects of finite cladding width on the dispersic I characteristics of optical fiber. Ackllllowledgements ‘l%e authors would like tc, thank Professor Cyrus C,antrell and Professor Wi liam Frensley for many 404.
15 June 1993
helpful discussions and support. This research was supported in part by U.S. Office of Naval Research grant NO0 14-92-J- 1030, and this support is gratefully acknowledged.
References [ 1] D. Gloge and E.A.J. Marcatili, Bell Syst. Tech. J. 52 ( 1973) 1563. [2] A. Gedeon, Optics Comm. 12 (1974) 329. [ 3 ] M. Matsuhara, J. Opt. Sot. Am. 63 ( 1973) 15 14. [4] A.W. Snyder, Electron. Lett. 6 (1970) 561. [ 51 T. Okoshi and K. Okamoto, IEEE Trans. Microwave Theory Tech. MMT-22 (1974) 938. [ 61 K. Oyameda and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-28 ( 1980) 1113. [ 71 SC. Gedam and S.S. Mitra, J. Lightwave Tech. LT-3 ( 1985) 706. [ 81 T. Tanaka and Suematsu, Trans. Inst. Electronics Comm. Eng.Jpn.E59(1976) 1. [9] C. Yeh and G. Lindgren, Appl. Optics 16 (1977) 483. [ lo] K. Okamoto and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-26 (1978) 109. [ 1 l] C.C. Su, Electron. Lett. 21 (1985) 858. [ 121 T.A. Lenahan, Bell Syst. Tech. J. 62 (1983) 2663. [ 131 J.D. Decontignie, 0. Parriaux and F. Gardiol. J. Opt. Comm. 3 (1982) 8. [ 141 C.C. Su and C.H. Shen, IEEE Trans. Microwave Theory Tech. MTT-34 (1986) 328. [ 151 T. Okoshi, Optical Fibers (Academic Press, New York, 1977). [ 161 A.W. Snyder, IEEE Trans. Microwave Theory Tech. MTT17 (1969) 1130. [ 171 D. Gloge, Appl. Optics 10 (1971) 2252. [ 181 G.D. Smith, Numerical solution of partial differential equations: finite difference methods (Clarendon Press, Oxford, 1979). [ 191 G. Strang, Linear algebra and its applications (Academic Press, New York, 1976). [ 201 J.H. Wilkinson, The algebraic eigenvalue problem (Oxford University Press, 1965). [21] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical recipes in C (Cambridge University Press, New York, 1988). [ 221 D. Gloge, Appl. Optics 10 ( 197 1) 2442. [ 231 K. Okamoto and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-24 (1976) 416. [ 241 T.I. Lukowski and F.P. Kapron, J. Opt. Sot. Am. 67 ( 1977) 1185. [ 251 E.K. Sharma, I.C. Goyal and A.K. Ghatak, IEEE J. Quantum Electron. QE-17 ( 1970) 23 17. [26] C.C. Su, IEEE J. Quantum Electron. QE-21 (1985) 1554. [ 271 S.I. Hosain, I.C. Goyal and A.K. Ghatak, Optics Comm. 47 (1973) 313.
99 ( 1993) 393-404
OPTICS
COMMUNICATIONS
Full length article
Analysis of optical fibers with arbitrary refractive index profiles: accuracy, convergence, and effects of finite cladding Lakshman
S. Tamil
Erik Jonsson School ofEngineering Richardson, TX 75080. USA
And Computer Science and
Centerfor Applied Optics, The University of Texas at Dallas,
and Gregory H. AickIen Erik Jonsson School of Engineering And Computer Science, The University of Texas at Dallas, Richardson, TX 75080, USA Received
24 September
1992; revised manuscript
received 4 January
1993
We have formulated a matrix eigenvalue problem for cylindrical optical fibers from a set of finite difference equations. Numerical solution of this problem yields the propagation constants for propagating modes. The method can be used for arbitrary index profiles, does not require the explicit evaluation of Bessel or modified Bessel functions, and does not use iterative methods to search for the propagation constants as was the case in earlier proposed methods using finite differences. The method is accurate, fast, and simple. We have established the convergence and stability of this method, and explored the effects of finite cladding width on the dispersion characteristics.
1. Introduction Wave propagation in optical fibers has been analyzed using various methods. We will be using a finite difference method. Other methods proposed to find the propagation constants of guided modes in optical fibers with arbitrary refractive index profiles include the WKBJ method, variational method, power series expansion method, staircase approximation method, and finite element method. The WKBJ method [ 1,2] is a geometrical optics approximation that works whenever the refractive index of fiber varies only slightly over distances of the order of the optical wavelength and are applicable only to thick fibers in which many modes can propagate. For those fibers in which only a few modes propagate, the error of the WKBJ method increases intolerably and this method is not applicable to modes near cutoff. Besides, the effect of an index valley at the core-cladding boundary, which plays an important role in reducing multimode dispersion, cannot be treated by the ordinary WKBJ method. In 0030-4018/93/%06.00
0 1993 Elsevier Science Publishers
the variational method the scalar wave equation is converted into a variational problem subject to the given boundary conditions. The variational problem is solved either by using the Rayleigh-Ritz method [ 31 or perturbation method [4]. In the RayleighRitz method the eigenfunction is expressed in terms of a set of orthogonal functions and the variational function is minimized. The disadvantage is that we need to assume a trial function [ 5 1. In the perturbation method of analysis, the computation of the propagation characteristics for an arbitrary profile is done by correcting the solution for a uniform core fiber considering the difference in the profile as the perturbation term. The power series expansion method [ 61 consists of expressing the refractive index for the field term by term. This method is useful only for cases in which the refractive index profile can be expressed by a relatively simple power series. In some cases the series do not converge and this method is not applicable [7]. In the staircase approximation method [8,9] the refractive index is approximated by an appro-
B.V. All rights reserved.
393
OPTICS
COMMUNICATIONS
staircase function. The ave equation is solved In each stratified layer and te solutions are then ‘.:onnccted at the cylindrical bt ndaries between these layers to obtain the proper sc ltion representing the propagation characteristics. he number of layers should be infinite in order tl .t the refractive index profile approaches that of act al fiber profiles. Thus, the results of the propagation instant will differ from the actual values when a fin e number of layers is used. A large number of laye requires considerable computer time and hence in .his method the accuracy and computer time are -aded off. The fiber problem has bee analyzed by Okamato and Okoshi [ lo] using a fin : element method formulated in the axial fields. ‘he problem with this method is that it suffers fror spurious modes when the finite elements are not arefully chosen [ 111. Lenahan [ 121 has formulat 1 a matrix eigenvalue problem from a finite elen :nt analysis using the Galerkin weighted residua method. To achieve computational efficiency, a 1 ecewise linear approximation to the solution func on must be used. In this paper we present 2 efficient finite difference method to find the proI gation constants of optical fibers with arbitrary rt ractive index profiles. The method does not involv a search procedure to find the propagation consta ts, nor does it require explicitly evaluating Besse and modified Bessel functions, as was the case in .he earlier works on finite difference analysis of 01 icalfibers [ 13,141. We construct a matrix equatior from a set of simultaneous finite difference ec .ations governing the propagation in an optical fi :r and solve for the eigenvalues to obtain the pr pagation constants. In sect. 2 we give the mathem; ical formulation of the discretized differential equat in at various grid points in the radial direction and t : construction of a matrix equation incorporating le boundary conditions at the core-cladding inteif :e and the jacket. Extending our method, which s formulated for a-index fibers, to arbitrary ref ctive index profiles is covered in sect. 3. In sect. 4 e discuss the numerical evaluation of propagation c nstants and present results. This includes a discu ion of the convergence and stability of the method ilong with the effect of the number of grid points I Lthe computation, and the effects of finite claddi g width on dispersion 1-n-k tr’
394
characteristics.
2. Mathematical
15 June 1993
Our conclusions
are given in sect. 5.
formulation
The optical fibers considered here are inhomogeneous dielectric cylinders of radius Q called the “core” surrounded by a homogeneous refractive index medium called the “cladding”. The cladding, in turn, is encased in a highly lossy material called the jacket. A representative fiber cross-section is shown in fig. 1. The refractive index profile of the fiber, called an a-index profile, is given by n2(1)=n:[1-2pd(r/a)“], =n:[1-2d])
forO
(1)
Here, n, is the maximum refractive index of the core, d is the relative refractive index difference between the core axis and cladding, p a parameter representing the refractive index step or valley at the corecladding boundary. A smooth continuation at the core-cladding boundary, the presence of a step, and that of a valley are expressed by p = 1, p < 1, and p > 1, respectively. {a 1CxeR} is a profile parameter. Some examples of I-u-index profiles are shown in fig. 2. The propagation characteristics of an optical fiber are governed by the scalar Helmholtz differential equation [ 15 ]
Jacket -L Fig. 1. Optical fiber showing grid points used in the example.
Volume 99, number
5,6
OPTICS
Core
15 June 1993
COMMUNICATIONS
Cladding
Region
1 .52
Step
rho Parabolic
Region
Index
1
=
Index
1.o Radius, r Fig. 2. Examples of a-index profiles. LY=co yields a step index, while (Y= 2 yields the parabolic index. Values of p control the characteristics of the interface between the core and the cladding; p< 1 results in a step at the interface, p> 1 yields a valley.
n2(r)k2-p2-
$
y=o.
This scalar wave equation is the simplification of the exact vector wave equation under the assumption that Vn/n is small, which includes the “small index gradient” and “weakly guiding approximations” [ 16,171. In the above equation v(r) is the transverse field function which may denote either the dielectric field or the magnetic field, r is the radial coordinate, n(r) is the radial refractive index profile, k is the vacuum wave number, fi is the propagation constant which is to be computed, and m is a mode parameter given by m=l,
forTEandTMmodes
(n=O),
=n+l,
for EH modes (n~n\l) ,
=n-1,
for HE modes (noN)
.
=O,
(2) Y(O)=0
We need to solve the differential equation in order to compute the propagation constants. From the rotational properties of w the associated boundary condition at the center of the core (r=O) is
formp0.
9
(4)
The other boundary condition applied tinction of field at the jacket written as vjacket
=
vru,= b =
0
is the ex-
>
(5)
where b is the radius of core and cladding
together.
2.1. Transformation to nondimensional form We need to nondimensionalize the differential equation for easy computation. This is achieved by setting u=wlv0,
(3)
form=O,
x=rla,
(6)
where v/Ois the maximum field amplitude and a is the radius of the core. Substituting eq. (6) into eq. (2) we obtain kzn2(xa)-/12-$$
>
u=o.
(7)
395
15 June 1993
OPTICS COMMUNICATIONS
by including the refractive i dex distribution given ly cq. ( I), the above equatil 1 can be rewritten as
( 11 ), and defining
-!Fy
CI’U
+a2 k2nT[ l-_/j
---; + C-Y-
h is the width between the grid points and x=ih, {i=O, 1, 2, . ..I. Substitutingeqs. (13) and (14) into
:a)]-/32-
?J
g
(
>
V2nT
17,
=: 0,
(15)
n:--nz’
(8)
we get
where
J‘(r)=2pA(r/a)U, =2A,
a
Defining
the parameters
,li=a(k2n:
-P2)
we can define
u~+l-2ut+ui_I + fUi+l-Ul_I
OQr
“2 )
h2
ed frequency,
propagatic
(9) I constant,
U-i=0 3
and on rearranging, as [ 171
V2= U2+ W2=k2a2(nT -ni and the modified
2h
+
C and Was
WE& P2-k2n:)‘/2,
V, the normali
rh
8, as
the equation
(16)
becomes
u,_i[-$(l-t>] +u,[$(2+ $) +Bliiha)-ij]
=a2(k2n:-/3’). +.,+,[-;(1+;)]=0.
(10) ‘Then eq. (8) becomes
p---V2n:
!e+;g+
(
?Z: -n:J
2 ax,->
u=o, 1
(11)
with
O,
f(x)=2pdxa,
(12)
x>l.
=2A,
2.2. Discretizing the differer; ial equation When the function u and is derivative are single functions of x, the first valued, finite and continuou an be approximated by and the second differentials third order difference form. las as follows [ 18 ] :
du
-z
2h
dx
d2u s=
(13)
’
h2
396
ui+,
+I?f(lha)-B)
+u&$
=O,
(18) where
=O,
m=O, m#O.
(14)
.
u* )
4+2m2-6 2h2
Ul
At i=2,
Ui+i -22ui+U,-i
Here ui=u(x)
For the purpose of illustration, we have chosen six grid points along the radial direction as shown in fig. 1. In general, the number of grid points can be any number not less than four, the minimum necessary to take care of the boundary conditions. Depending on whether m = 0 or m # 0, the field or its derivative vanishes at the center of the core. When the derivative of the field vanishes, u0 = u,. Writing finite difference equations at the grid points, we obtain the following set of equations. At i= 1,
6=1,
%+I -u,-I
(17)
=u(x+h)
u,_, =u(x-h)
,
(3 +u,($ =o.
+Ff(2ha)-j?)
+u&$ (19)
Volume 99, number
$6
OPTICS
COMMUNICATIONS
At i=3,
15 June 1993
- (2i+ 1) ai,,+
+Yf(3ha4
+u,($)
=o.
(20)
1 =
2ih2
(25)
’
In order to convert the problem into an eigenvalue problem, we rewrite eq. (24) as (T-fl)u=O,
At i=4,
=o.
(21)
At i=5, since the field goes to zero at the jacket,
(22) Finally, at i= 6, again using the boundary that the field goes to zero at the jacket,
(26)
where I is the identity matrix, and T is a tri-diagonal matrix. Equation (26 ) defines an eigenvalue problem. This means that eq. (26) has a nontrivial solution if and only if p”is an eigenvalue [ 191. Hence, the required normalized propagation constants contained in Fare obtained by finding the eigenvalues of the tri-diagonal matrix T. This mathematical formulation can be generalized to {n 1n 1EN} grid points in the radial direction of the fiber without difficulty.
condition
3. Arbitrary profiles, multipie layers, and field ug=o.
(23)
Since the boundary condition in eq. (23) is incorporated in eq. (22), we have a system of five equations. 2.3. Matrix equation formulation Formulating a matrix equation from the above set of equations for the convenience of generalization and easy computation, we obtain a,, -s a21
a12 a22
Au=
Ul
-P
a32
a23 a33
u2
-P
a43
a34 a44
-
u3 B
a54
i
a45 a55-B
=o,
iii
~5
(24)
where the matrix elements
ai+, = a,,i =
u4
- (2i-
are defined
by
1)
2ih2
2i2-!-2m2-6 2i2h2
’
+pf(iha),
2i2+m2 = i~h2 + pf(iha)
,
i=l
i#l,
,
distributions We have developed this method of analyzing optical fibers using the a-index profile. This is because the a-index profile is commonly used in the literature and can represent a large variety of real refractive index profiles, including the very important step and parabolic profiles. But our formulation is not limited to a-index profiles. To see how to extend the method to arbitrary profiles without rederiving a system of finite difference equations, consider eq. ( 11). The refractive index profile is included in this equation through the functionf(x), which is defined in eq. (12). Usingf(x), the refractive index profile of the fiber can be rewritten as
n2(x)=n:[
1 -f(x)]
Solving for f(x) f(x)=
1 -n2(x)/nf.
.
(27)
yields (28)
By generating the discretized A=f(x,= ih ) from samples of an arbitrary refractive index profile n(x,), the method we have outlined in this paper can be used directly on arbitrary profiles, as long as the “weakly guiding” approximation holds. Multiple layer waveguides of any number of layers may be considered special cases of arbitrary refrac-
397
‘~OlUllli 0
1’.
nurlbcr 5.‘)
OPTICS COMMUNICATIONS
t ive inde!x profiles Since WC’II Lve normalized the fiber c’#:)rcradius. ,z, to unity, ‘rv~’must explicitly define which layer:, comprise the )re before using our methcd and scale all quantit s accordingly. Since the propagation const mt for a mode i, /I$, is uniquely associated with a fi .he field distributions for a propag.sting mode can k determined from the cigenve’ctors u of matrix T ( I( e eq. (26) ). Many eigenvalue routines will return (’ genvectors as well, but at the cost of greatly incr’: lsing the number of computations. Useful approximations tcl the eigenvectors for propagating modes can be co. lputed by constructing the tri-diagonal matrix T, SI btracting a specific fl from each element of the ma:’ Ldiagonal, and solving for the melem:nts of u using st;l ldard techniques from linear algebra. From the ob:l rvation that for propagating modes the field will approach zero at the cladding/jacket boundary, V’ can set uN, the rightmost element of u, to a very ! mall value (not zero), and use a simple backsubst ution process to solve for the rest of the u,. This 11ocedure yields a good approximatj on to the field di! tribution multiplied by an arbi trar), constant.
4. Nwneriad
evaluation,
rest Its, and discussion
Although the mathemati: al formulation of our method for determining the propagation characteristics ad an optical fiber is co Lched in terms of matrix equations, there are special structures that lead to very efficient numerical imp mentations. First, since T is a tri-diagonal matrix, vi : can use sparse matrix techni’ques to reduce storal e requirements for T. Second, since T is a quas symmetric tri-diagonal matrix, we can use a sirnil lrity transformation to convert T into a real, symt metric matrix [20]. Finally, ,the ei genvalues of a re: I , symmetric matrix may be computed using an efficil: It 0 (N2 ) algorithm (in our case, the tq1i.c routine f [ ,rn ref. [ 2 11, which has an alperation count of apprr ximately 30N’). Using eqs. (25), we havt implemented a pair of C language: programs which :ompute the normalized propagatic’n constants for ! ibers with arbitrary refractive index profiles. WI: define the normalized propagation constant as 398
x_
;z _
15 June 1993
k2n?-D2 k*n: -k*n:
’
Note that some authors normalized propagation b=1-U2/V2=1-x.
(29) (e.g. Gloge [ 221) define a constant as (30)
One program computes x for all propagating modes at a specific value of m in eq. (2) over a range of normalized frequencies V. Another program searches for the cutoff frequency ( V,) of a specific linearly polarized (LP) mode. Both programs allow all of the parameters in eq. ( 1) to be varied, as well as the values of b and N,, which are the fiber radius (see fig. 1) and number of grid points in the fiber core, respectively. In verifying the performance of our method, we have computed the propagation characteristics of step index and parabolic index fibers over a normalized frequency range of 0 to 20. These index profiles have analytical solutions and have been studied analytically and numerically by other authors [ 14,23-261. Our results agree well with previously published results, as shown in table 1. Note that for propagating modes, x must lie between 0 and 1 (i.e. O
Volume 99, number
56
OPTICS
COMMUNICATIONS
15 June 1993
Table I Comparison of the cutoff frequencies obtained by the finite difference method with analytical and previous numerical results. Raber is the fiber radius, R,, is the core radius, and 6 is the percentage difference from the infinite cladding result. 256 points were used in the fiber core. cr
1 2
Infinite cladding
Mode (m, I)
4.381 3.518 7.45 1 5.744 9.645 7.848 9.904 3.181 3.000 2.886 2.649 2.527 2.405
1,l I,1 12
2,1 22
3 4 5 10 20 co
Normalized
3,l 4,l 1,l 1,l 171 1,l 1,l 121
1
cutoff frequency
R m-x= lOR,re
6 (%)
Rfiber
4.391 3.526 7.457 5.744 9.645 7.848 9.904 3.189 3.007 2.894 2.657 2.535 2.413
0.23 0.23 0.08 c lo-* < 1o-2 <10-z < 1o-2 0.3 0.2 0.28 0.30 0.32 0.33
4.384 3.520 7.453 5.744 9.645 7.848 9.904 3.183 3.002 2.888 2.65 1 2.529 2.407
=
2’%0,
6 (%) 0.07 0.06 0.03
.o
5
Fig. 3. Dispersion
characteristics
15
10
Normalized
Frequency,
of a step index fiber (a=a,
frequency using from 4 to 256 points in the core, and for fiber radii from 5 to 20 times the core radius. From this figure we can see the expected convergence on a final result as the number of points in the
20
V
A=O.O38). R,+._=IOR,,.
core increases. The effect of cladding width is also apparent. The effect of the number of grid points in the core is two fold. As the number of grid points is increased the distance between samples of the refractive index
399
‘L’ol~.lmc9’7, number 5,b
OPTICS
COMMUNICATIONS
5
10
Normalized Fig. 4. Dispersior
0.0
characteristics
c
’
’
’
I
I
I
I
I
5
400
characteristics
V
I,
I
I
10
Normalized Fig. 5. Dispersion
of a parabolic
20
15
Frequency,
index fiber (cu=2, p= 1,4=0.038).
ofa parabolic
’
15 June 1993
Frequency,
I
I 15
I
I
Rfikr= IOR,,.
!
I
I 20
V
index fiber (cy=2, p=O.75, d=O.O38). RBkr= IOR,,.
Volume 99, number
5,6
OPTICS
COMMUNICATIONS
15 June 1993
L
0.0
0
5
10
Normalized Fig. 6. Dispersion
characteristics
’ ’
-0.05 0
’
’
’ 50
ofa parabolic
’
’ ’ ’
’
100
’
15
Frequency,
20
V
index fiber ((r=2,p=2,4=0.038).
’ ’ ’
’ 150
’
’ ’ ’
’ 200
’ ’ ’
Rfibcr= lOR,.
’
’ 250
’
’
’ ’ 300
No. of Grid Potnts in the Core Fig. 7. Convergence behavior of computed cutoff frequency with the number of grid points in the fiber core for the LP,, and LPO* modes of a step index fiber. 401
Volume 99. number 5,6
OPTICS COMMUNICATIONS
profile is decreased, resulting in a better approximation of the actual profile. This is especially apparent in profiles with sharp t ansitions at the corecladding interface, such as fc r a!>> 1. Cases where p # 1 (see eq. ( 1) ) are also lil ely to be poorly modeled by a small number of co *e sample points. The effect of reducing the number if sample points in the core is to apply a “low pass” spatial filter to the refractive index profile. Setting the number of poir ts in the core also effectively applies a filter to the spatial frequency content of the field distributior s calculated for each mode in the fiber. When ccImputing propagation constants at higher norrnalizc d frequencies, using a small number of samples ma t induce errors due to a form of “aliasing”. These 1wo effects are responsible for the poor results wh :n the number of grid points in the core is below ap proximately 10 for the modes we have examined. Using the step index fiber as an example, fig. 8 demonstrates the effects of th : number of grid points by plotting the computed cute ff frequency for modes LPO,, LP, ,, and LP,, for sew ral different grid sizes. In this figure, each curve is r ormalized to the value
100
15 June 1993
of the cutoff frequency for that mode calculated with 256 points in the core. We can see that for mode LPol the cutoff frequency calculated with 8 points in the core is less than 1.0015 times that computed using 256 points in the core, while for mode LPI, (with a higher cutoff frequency) we need at least 12 points in the core for similar results. In general, as the normalized frequency increases, the number of points in the core must be increased to maintain the accuracy of the method. For modes with relatively low cutoff frequencies, variations in cladding width produce large changes in the calculated cutoff frequency, V,. Cutoff frequency increases as the cladding width decreases. This is the expected behavior. Analyses assuming infinite cladding width, while adequate for many purposes, fail to account for the increasing importance of finite cladding width as the normalized frequency becomes smaller. The fundamental mode, which has no cutoff frequency when the cladding is infinite, shows a definite cutoff in real fiber. Figure 9 shows the effects of cladding width on the cutoff frequencies of two LP modes in a step index fiber. In this plot, the curves for each mode are nor-
150
200
No. of Grid Points in the Core
Fig. 8. Effect of the number of gri i points on the computed cutoff frequencies of propagating modes of a step index fiber (LPO,, LPI,, and LP,,) .
402
15 June 1993
OPTICS COMMUNICATIONS
Volume 99, number 5,6
Cladding Width/Core
Width
Fig. 9. Effect of cladding width on cutoff frequency for a step index fiber.
0.0
0.5
1.5
1 .o
Radius,
2.0
I
Fig. 10. Field patterns for the LP 01, LP,I, and LPtl modes of a step index fiber with 4~0.038. These patterns were computed using R tibcr=1ORc.m V= 10,and 32 points in the core.
‘I!) lr:, : j’i. number
5.6
OPTICS
COMMUNICATIONS
#‘II:1):’c 8.110 the cutoff frequent for that mode in the nl 1.1‘1t .: cladding case. As an ex lmple, consider mode 1.J , , 14/t:can see that the cum ‘ffrequency when the I‘ilzer T.rdius is 10 times the COI : radius is only about I .(rO ; imes (0.3%) greater thz t the cutoff frequency wt. er I.hazfiber radius is 25 t mes the core radius. How e #Tel-,when the fiber radi ts is only 5 times the core mdius then the cutoff f equency increases to 1.O14 I.imes ( 1.4%) greater th: n the cutoff frequency when Ime fiber radius is 25 ti les the core radius. Usi’rg the approach outline d in sect. 3, we have camp .ned the field distributic IS for three LP modes in a !il repindex fiber. Figure 10 shows the results, with the -fis:ld patterns normalized so that the maximum VElue in each pattern is one. 1 hese results agree well u’ th, I he results reported in r-1f. [ 15 1.
5. Co ,dlusions We know that when d is sn all in an optical fiber, tit e s( al;ar approximation yiel Is results that are very close ‘to the exact vector forrr ulation. Even for large dlIfft:r lences between the core .nd cladding refractive itIde> es, optimum single-mot e fiber parameters obt;.inel:I ffrom the scalar apprc ximation differ negligjbly from those obtained u! _ng the exact formulation [ 271. V& heave developed a methl d to evaluate the propagatii:m constants by transfc .ming the scalar wave equa ;ion into a set of finite d fference equations and then converting into a matr x eigenvalue problem. The imethod does not involv a search procedure to l‘nd I.he propagation constat ts, or the explicit evalIIatic n of Bessel and modified Bessel functions, which i:s tirae consuming, as was tf : case in earlier works. W;: lhave demonstrated t te convergence of the tnetl~od and the dependence of the rate of conver,g;ence on the number of grit points. The method is eleg; nt., stable, straight forw .trd, is applicable to arbitrary index profiles, and ir accurate. We have also expllored and established the effects of finite cladding width on the dispersic I characteristics of optical fiber. Ackllllowledgements ‘l%e authors would like tc, thank Professor Cyrus C,antrell and Professor Wi liam Frensley for many 404.
15 June 1993
helpful discussions and support. This research was supported in part by U.S. Office of Naval Research grant NO0 14-92-J- 1030, and this support is gratefully acknowledged.
References [ 1] D. Gloge and E.A.J. Marcatili, Bell Syst. Tech. J. 52 ( 1973) 1563. [2] A. Gedeon, Optics Comm. 12 (1974) 329. [ 3 ] M. Matsuhara, J. Opt. Sot. Am. 63 ( 1973) 15 14. [4] A.W. Snyder, Electron. Lett. 6 (1970) 561. [ 51 T. Okoshi and K. Okamoto, IEEE Trans. Microwave Theory Tech. MMT-22 (1974) 938. [ 61 K. Oyameda and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-28 ( 1980) 1113. [ 71 SC. Gedam and S.S. Mitra, J. Lightwave Tech. LT-3 ( 1985) 706. [ 81 T. Tanaka and Suematsu, Trans. Inst. Electronics Comm. Eng.Jpn.E59(1976) 1. [9] C. Yeh and G. Lindgren, Appl. Optics 16 (1977) 483. [ lo] K. Okamoto and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-26 (1978) 109. [ 1 l] C.C. Su, Electron. Lett. 21 (1985) 858. [ 121 T.A. Lenahan, Bell Syst. Tech. J. 62 (1983) 2663. [ 131 J.D. Decontignie, 0. Parriaux and F. Gardiol. J. Opt. Comm. 3 (1982) 8. [ 141 C.C. Su and C.H. Shen, IEEE Trans. Microwave Theory Tech. MTT-34 (1986) 328. [ 151 T. Okoshi, Optical Fibers (Academic Press, New York, 1977). [ 161 A.W. Snyder, IEEE Trans. Microwave Theory Tech. MTT17 (1969) 1130. [ 171 D. Gloge, Appl. Optics 10 (1971) 2252. [ 181 G.D. Smith, Numerical solution of partial differential equations: finite difference methods (Clarendon Press, Oxford, 1979). [ 191 G. Strang, Linear algebra and its applications (Academic Press, New York, 1976). [ 201 J.H. Wilkinson, The algebraic eigenvalue problem (Oxford University Press, 1965). [21] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical recipes in C (Cambridge University Press, New York, 1988). [ 221 D. Gloge, Appl. Optics 10 ( 197 1) 2442. [ 231 K. Okamoto and T. Okoshi, IEEE Trans. Microwave Theory Tech. MTT-24 (1976) 416. [ 241 T.I. Lukowski and F.P. Kapron, J. Opt. Sot. Am. 67 ( 1977) 1185. [ 251 E.K. Sharma, I.C. Goyal and A.K. Ghatak, IEEE J. Quantum Electron. QE-17 ( 1970) 23 17. [26] C.C. Su, IEEE J. Quantum Electron. QE-21 (1985) 1554. [ 271 S.I. Hosain, I.C. Goyal and A.K. Ghatak, Optics Comm. 47 (1973) 313.