Variational and finite element analyses of nonlinear strip optical waveguides

Optics Communications 94 (1992) 37-43 North-Holland

OPTICS COMMUNICATIONS

Variational and finite element analyses of nonlinear strip optical waveguides Q.Y. Li ~, R.A. S a m m u t a n d C. P a s k Department of Mathematics, University College, University of New South Wales, Canberra, Australia Received 2 January 1992; revised manuscript received 17 July 1992

A new, iterative finite element scheme for solving nonlinear wave guiding problems is proposed. Unlike previous numerical methods, this approach gives both the stable and unstable stationary solutions of a nonlinear strip waveguide. We use a double gaussian trial function to solve the same problem by variational principles. This approximation gives very good results and provides a simple model to characterize nonlinear stationary waves. The influence of waveguide parameter on nonlinear behaviour is discussed. We also apply both the variational and finite methods to analyze the stripe surface waves at a linear-nonlinear interface.

1. Introduction N o n l i n e a r optical waveguides offer a wide range o f potential applications, one o f the most i m p o r t a n t being all-optical switching. F o r example, optical switching has been p r e d i c t e d [ 1 ] in a slab wavequide with a linear core a n d n o n l i n e a r cladding. At low power, the field in this waveguide is confined m a i n l y to the core b u t as p o w e r increases, the field m a x i m u m moves out o f the core a n d into the nonlinear claddding. In the past, most attention has been c o n c e n t r a t e d on o n e - d i m e n s i o n a l solutions o f slab waveguides. Recently, however, spatial solitons have been experimentally observed in slab waveguides with a nonlinear core a n d linear cladding [2]. These are fields in which c o n f i n e m e n t is p r o v i d e d in one transverse d i m e n s i o n by the linear refractive index difference a n d in the orthogonal d i m e n s i o n b y the nonlinear self-trapping o f the beam. Both n u m e r i c a l [ 3 ] a n d simple, v a r i a t i o n a l solutions [4] have been presented for the steady-state b e h a v i o u r o f these solutions. H o w e v e r these waves are not confined in b o t h transverse d i m e n s i o n s at low powers so, for practical

On leave from the Northern Jiao Tong University, Beijing, China.

applications, there m a y be advantages in using a twod i m e n s i o n a l guiding structure. Chrostowski a n d Chelkowski [ 5 ] a n d Ettinger et al. [ 6,7 ] discuss structures consisting o f linear rib or strip waveguides with either a nonlinear substrate or nonlinear cover. These p r o v i d e guidance in b o t h transverse directions and m a y be practical configurations for optical switching devices. To solve for waves in relatively c o m p l i c a t e d structures like this, Ettinger et al. [6,7] use an iterative finite element m e t h o d in which a sequence o f linear p r o b l e m s is solved with the solution at each step being used to generate the refractive index profile for the next step. Their iteration scheme, like the one used in one-dimensional finite element analysis [ 8,9 ], fails to produce solutions within and near the unstable region. This is interpreted as illustrating a link between the numerical process for building up the i n d u c e d index profile a n d the underlying physical process. In this paper, we use a commercially available finite element package a n d a different iteration scheme to analyze a waveguide c o m p o s e d o f a linear strip o f height a and width b sitting on a nonlinear substrate a n d s u r r o u n d e d by a linear cladding as shown in an inset in fig 1. By our approach, b o t h stable a n d unstable solutions can be found a n d we are able to accurately predict the two turning points which are the

0030-4018/92/$05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.

37

Volume 94, number 1,2,3

o

OPTICS COMMUNICATIONS

1 November 1992

[k2n2(x, y) +k2f(au] 2) -f12]~,=O

02__~ ~ + a2_~_ ~ + Ox2 Oy2

a

(1)

tlcl

¢.o

b = 0 ' 4 , ,''/''

t.t)

i

I

0.0

I

I

0.04

I

I

0.08

I

0.12

I

0.16

where ~,(x, y) is the scalar field, fl the propagation constant, k the wave number n(x, y) the linear refractive index distribution and the function f(a~/2 ) describes the nonlinearity, a=n~n2ceo where ns is the refractive index of the substrate and n2 is the nonlinear Kerr coefficient of the medium. Initially, we consider a nonlinear material with gaussian saturation characterized by

POWER [ rnW ]

Fig. 1. Effective index as a function of power by the infinite element method, showing the influence of the strip height b. Strip width, a=2t~m, b= 1.2 vtm (solid line) and b=0.4 vtm (dashed line). Other waveguide parameters are given in the text. The inset figures illustrate the waveguide structure and typical field distribution.

2GAnsat/A ' where Ansat is the maximum nonlinear index increase. For weak fields, this gives a Kerrqike nonlinear behaviour. Now let

~(x, y) =a¢(x, y), most important parameters in characterizing the switching behaviour. In waveguide analysis and design, it is always desireable to have simple approximate solutions which characterize the main features of the problem under consideration. This is particularly true for the problems discussed in this paper, because strict numerical analysis is relatively difficult. For example, Chrostowski and Chelkowski [ 5 ] use a modified effective-index method. We have successfully used a one-dimensional variational method to produce accurate solutions and even to predict stability behaviour of nonlinear waves [ 10,11 ]. A two-dimensional variation method using such simple trial functions as gaussian and sech functions also proved to be successful in describing spatial solitons in a slab waveguide with nonlinear core [4 ]. In this paper we use a three-parameter double gaussian trial function in a variational expression and again obtain very good results.

2. Finite element analysis In this section, we find the stationary solutions of the nonlinear, scalar wave equation

38

where

O(x, y)

is a normalized field satisfying

i i (~2(x,y) dxdy=l. --oo

(3)

(4)

-oo

The modal power is then given by

p= tic%A2

(5)

2k'and the function ¢ satisfies 02~ + 02___~ ~ + OX2

O.y2

[k2n2(x,y)-fl 2

+k2f(aA2¢2) 1•=0.

(6)

In the existing iterative, finite element scheme [6,7 ], an equation equivalent to eq. (6) is solved at a specific power by assuming initially that the field in the nonlinear term is that of the linear mode. A finite element method is then used to solve for the mode of the linear waveguide with refractive index profile n 2+ f ( a A 2¢~2). Iterating this process eventually leads to a self-consistent solution for the nonlinear problem. However the method only converges for the stable modes. We propose a modification to this method in which we assume a given fl and then solve for the appro-

Volume 94, n u m b e r 1,2,3

OPTICS C O M M U N I C A T I O N S

1 November 1992

priate A and ~ "1. Our iteration scheme is shown in the following equation, where we regard the field magnitude, Am+ 1, as the eigenvalue. 02~m+l .~ 0 2 ~ m + l OX2 Oy 2 "1-4 2 + 1

( k2f(

+Fk2rl2(x,y)_fl2 l 2° x7 )J~m+l

=0

.

(a)

(b)

(7)

Starting with an initial guess of ~o and Ao, we solve eq. (7) iteratively using a commercially available finite element package [ 13 ]. At the end of each step, is normalized according to eq. (4). We assume symmetry in the x direction to reduce the size of the problem and impose a zero boundary condition far from the field center (typically at a distance 2 ½-3 times the dimensions of the guide ). About 1000 linear elements are used with larger elements near the boundary and the mesh moving with the field peak to improve efficiency. We have checked this iteration scheme against the analytical solutions for one-dimensional problems and found excellent agreement. In all the numerical examples in sects. 2 and 3, we assume the values nco= 1.57, nd= 1.55, ns= 1.55, wavelength 2=0.515 ~tm, n2= 10 - 9 m 2 / W (MBBA liquid crystal) and, except in fig. 4, Ansa t = 0. l. Figure 1 shows numerical results for a fixed strip width ( a = 2 ~tm) and two different strip heights. For b = 0.4 ~tm (dashed line), the field gradually moves out of the core as power increases, while for b = 1.2 ~tm (solid line), there is a quick switching and bistability appears. Our results give both stable (positive P-fl slope) and unstable (negative P-fl slope) solutions without any ambiguities. Figure 2 shows the field distribution contours corresponding to points A and B in fig. 1. As we would expect, the field on the lower branch is approximately the linear solution and is mainly confined in the core, while the field on the upper branch moves to the nonlinear substrate. Figure 3 illustrates the influence of the strip width, a, on the solution. As expected, switching occurs at a lower power in the case of the narrower strip width because the field is less strongly bound to the core. ~l A similar procedure has recently been applied to a vector finite element analysis by Wang and Cambrell [ 12 ].

Fig. 2. Field contours calculated using the finite element method. (a) The field corresponding to point A in fig. 1 where the field is mainly confined in the core and (b) the field corresponding to point B where the field has moved to the nonlinear substrate.

bI

ncl ns

-x: CO LO

,,/~

~

S/

,.2 . . . . . . . . . . . . . .

0.0

......

a = 2pro =

"

0.04

0.08

O. 12

0.16

POWER [ mW ] Fig. 3. Effective index as a function of power by the finite element method, showing the influence of the strip width a. Strip height b = 1.2 ~tm, a = 2 v.m (solid line) and a = l ~tm (dashed line). For smaller a, the field is more concentrated and switching occurs at a lower power.

o

Ansat = -~ CO

~u-t.

0.21// //no \\

_0 02

....~

Ansat = 0.02

tD

0.0

0.04

0.08

0.12

0.16

POWER [ mW ] Fig. 4. Effective index as a function of power by the finite element method, showing influence of nonlinear saturation, a = 2 ~tm and b = 1.2 Ixm.

39

Volume 94, n u m b e r 1,2,3

OPTICS C O M M U N I C A T I O N S

Finally we look at the influence of saturation in fig. 4. As we would expect, when the nonlinearity saturates at low levels, the behaviour of the waveguide is very nearly linear whereas when the maxim u m nonlinear index change is sufficiently large, a self-trapped wave exists above a certain threshold power.

3. Variational analysis

The nonlinear scalar wave equation can be transformed into the variational problem of finding a stationary point of the functional [ 10 ]

1 November 1992

sen so that constraint (10) is satisfied. By substituting the trial function (11 ) into the variational expression (8), and optimizing J with respect to the three parameters tox, toy and Yo, we can find the best possible approximation. Figure 5 shows the P-fl relationship by both the finite element method (solid line) and the variational method (dashed line) for the two strip heights used in fig. 1. The variational method gives accurate solutions in the stable region, while in the unstable region, attempting to optimize J does not give any results because, in that region, the solution of eq. (1) corresponds to a saddle point rather than a maximum [ 11 ]. As we shall see in a later example, variational solutions can still be found in the unstable region

ii

J=

+ k2nZ(x, y)~,2+kEg(ot~,2)l dx dy

(8)

g

subject to the constraint

i i ~u2(x,y)dxdy=C. --oo

(9)

--oo

In eq. (8),

g

S

!

Fo.,e Var,a,,oo E. a

OOg2

z 1 2ns Ansat g(o~ )=~ f f(~)d~= o~ t°t~2-2nsAns"~

I

I

I

I

I

I

I

I

I

0

× [1-exp( =0,

°t~z ] ] ~ , 2n-~sat/_l)

for y < 0

fory>0,

(b)

(10)

,

~

and the constant C is related to power by

flC~oC~ ½nsc~oC. P= --~ We choose as our trial function a three-parameter double gaussian function defined as I

~(x, y) =

exp

2tox 2

2mr2 ] , (11)

where tox and toy account for the field spreading in x and y directions and Yo characterizes the movement of the field maximum. The amplitude is cho40

0.0

I

0.04

I

I

0.08

I

0.12

0.16

POWER[mW] Fig. 5. Effective index as a function of power by both finite element method (solid lines) and variational method (dashed lines).

Strip heights are (a) b=0.4 gm and (b) b= 1.2 gm. Strip width a=2 gm for both cases.

Volume 94, n u m b e r 1,2,3

OPTICS C O M M U N I C A T I O N S

1 November 1992

CJ

provided we look for stationary points rather than optima. This feature can be used to predict the stability of nonlinear stationary waves. (Incidentally, the arrows in fig. 5b indicate the behaviour of the variational solution when a step-by-step technique is used to calculate the curve and not necessarily the behaviour of the physical wave. We begin with the linear solution at low power and find the solution at subsequent power levels using the previous solution as a starting point. Once we have reached the turning point above 0.8 mW, a further slight increase in power leads to an optimum value on the upper branch of the curve. Similarly, if we reverse the process from high power levels, an abrupt jump in the variational solution occurs at the turning point near 0.5 mW. ) Next, we examine how the three variational parameters, Yo, cox and COy change with power. For b = 0 . 4 ~tm, which gives a monotonic P-fl curve, the three parameters also vary monotonically as shown by the dashed lines in fig. 6a. (In these figures, the solid lines give estimates of Yo, COxand COybased on fitting gaussians to the finite element solutions. ) The position of the field maximum, characterized by Yo gradually moves from the center of the strip to the nonlinear substrate. The field widths, COxand COy,first decrease quickly with power because of the nonlinear self-trapping effect and then increase slowly with power due to the nonlinear saturation effect. For b = 1.2 pm, where bistability appears in the P-fl curve, changes in Yo, COxand COywith power also show bistable behaviour as illustrated in fig. 6b. Once again, the arrows indicate the behaviour of the numerical solution.

0

c; 0

d

o,

9" i

i

i

i

o.

d

(t) x

~x cO

d

,5

d i

i

COy

0.0

i

i

0.05

0.10

0.15

POWER [mW ] (a)

0.0

i

i

0.05

0.10

0.15

POWER [ mW ] (b)

Fig. 6. Variational parameters Yo, tox and toy as functions of power with strip width a = 2 ~tm. Strip heights are (a) b = 0 . 4 ~tm and (b) b = 1.2 ~tm. Dashed lines give the parameters found from the variational method while solid lines are found by fitting a gaussian to the finite element solution.

the refractive index profile is defined by n2= no2+f(ot~,2) , y < 0 ,

4. Stripe surface waves Nonlinear waves of a very similar form to the guided waves discussed above have recently been shown [ 14 ] to propagate at the interface between two semi-infinite media - one linear and the other nonlinear. However, the numerical approach used in ref. [ 14 ] is an iterative FFT-finite difference procedure which is rather difficult to implement. In this section, we show that the iterative finite element method and the variational approximation can both be used to investigate these waves also. We use the same parameters as in ref. [ 14 ] where

=n~,

y>0,

(12)

and, in this case, the nonlinearity is produced by a two-level saturable system with ag/2 f(°~v/2) = 1 + / t a ~ 2 '

(13)

where/z = 1/ (2noAn~t). The linear refractive indices of the two media are nL=1.635 and no= 1.627, respectively. The finite element solution for this problem proceeds precisely as for the strip waveguide. The solid lines in fig. 7 show the calculated effective indices as 41

Volume 94, number 1,2,3

OPTICS COMMUNICATIONS

t..O

}.t= 0.3

£o £o

/, //

,,_.:

5. Conclusion

I I

I I I [ .//// "~

10 ~= 3.0 /

I

I

4

6

8

Power [roW] Fig. 7. Effective index as a function of power for the stripe surface wave at a linear-nonlinear interface. Parameters are the same as in ref. [ 12]. Solid lines are calculated using the finite element method and dashed lines using the variational approximation. functions o f power for three different saturation levels. In principle, the variational solution also follows the same p r o c e d u r e as above, but now the function g in eqs. ( 8 ) a n d ( 1 0 ) must be replaced by [ 10 ] g(o~/2) = ~ =0,

(o~q/2-- 1 in ( 1 + # o t ~ 2 ) ) , y>0,

y<0,

In this paper, we have p r o p o s e d a new iteration scheme for finite element analysis o f nonlinear waveguides. This scheme can find both stable and unstable stationary solutions, enabling us to accurately predict the two turning points on the bistability curve. The m e t h o d can be used to analyze nonlinear stationary waves in waveguides with arbitrary structure and arbitrary nonlinearity. Using this technique, we have investigated the influence o f waveguide d i m e n s i o n s a n d saturation on the b e h a v i o u r o f a nonlinear strip waveguide. We have also used a double gaussian trial function and solved the same p r o b l e m using a variational principle to obtain a simple, accurate a p p r o x i m a t i o n . Finally, we have shown that both the finite element m e t h o d and the variational a p p r o x i m a t i o n can be used to analyze stripe surface waves at a l i n e a r - n o n l i n e a r interface.

Acknowledgements We are grateful to the Australian Research Council and Telecom Australia for financial support.

(14)

Although they do not p e r f o r m a detailed stability analysis, A k h m e d i e v et al. [14] suggest that, by analogy with t w o - d i m e n s i o n a l nonlinear surface waves, these stripe waves are stable along the positively-sloping branches o f the P - f l curves a n d unstable along the negatively-sloping branches. O u r variational calculations support this suggestion. The function J has a m a x i m u m value along the positively-sloping branches only. However, as we found in the case o f the p l a n a r waveguide [ 11 ], it is possible to find solutions to the v a r i a t i o n a l p r o b l e m along the whole o f the P - f l curve if we look for stat i o n a r y points o f J rather than maxima. The d a s h e d lines in fig. 7 show the results. As can be seen, the variational m e t h o d is good enough to m a k e some qualitative predictions in this case a n d could be imp r o v e d by using a slightly m o r e complex trial function than the gaussian in eq. ( 11 ).

42

1 November 1992

References [1] See for example, G.I. Stegeman, C.T. Seaton, W.M. Hetherington, A.D. Boardman and P. Egan, in: Nonlinear optics: materials and devices, eds. C. Flytzanis and J.L. Oudar (Springer, 1986) p. 31; D. Mihalache, M. Bertolotti and C. Sibilia, in:Progress in optics, Vol. 27 (Elsevier, 1989) p. 229; A.D. Boardman, P. Egan, F. Lederer, U. Langbein and D. Mihalache, in: Nonlinear surface electromagnetic phenomena, eds. H.-E. Ponath and G.I. Stegeman (Elsevier, 1991) p. 73. [2] S. Maneuf, R. Desailly and C. Froehly, Optics Comm. 65 (1988) 193; J.S. Aitchison, A.M. Weiner, Y. Silberberg,M.K. Oliver, J.L. Jackel, D.E. Leaird, E.M. Vogel and P.W.E. Smith, Optics Lett. 15 (1990)471. [3] N.N. Akhmediev R.F. Nabiev and Yu.M. Popov, Optics Comm. 69 (1989) 247. [4] Q.Y. Li, C. Pask and R.A. Sammut, Optics Lett. 16 ( 1991 ) 1083; R.A. Sammut, C. Pask and Q.Y. Li, in:Technical digest on Nonlinear guided-wave phenomena 1991 (Optical Society of America, Washington, 1991 ) vol. 15, p. 399.

Volume 94, number 1,2,3

OPTICS COMMUNICATIONS

[5] J. Chrostowski and S. Chelkowski, Optics Lett, 12 (1987) 528. [6] R.D. Ettinger, F.A. Fernandez, B.M.A. Rahman and J.B. Davies, IEEE Photon. Technol. Lett. 3 ( 1991 ) 147. [7] B.M.A. Rahman, F.A. Fernandez, R.D. Ettinger and J.B. Davies, in: Technical digest on Nonlinear guided-wave phenomena 1991 (Optical Society of America, Washington, 1991 ) Vol. 15, p. 96. [8] B.M.A. Rahman and J.B. Davies, Int. J. Optoelectronics 4 (1989) 153. [9]K. Hayata, M. Nagai and M. Koshiba, IEEE Trans. Microwave Theory Tech. 36 (1988) 1207.

1 November 1992

[10] R.A. Sammut and C. Pask, J. Opt. Soc. Am. B 8 (1991) 395. [ 11 ] R.A. Sammut, Q.Y. Li and C. Pask, J. Opt. Soc. Am. B 9 (1992) 884. [12]X.H. Wang and G.K. Cambrell, Proc. 16th Australian Conference on Optical fibre technology, ( 1991 ) p. 314. [ 13 ] Finite Element Library, Rutherford Appleton Laboratory, Oxfordshire, UK (1982). [ 14] N.N. Akhmediev, R.F. Nabiev and Yu.M. Popov, Solid State Comm. 66 (1988) 981.

43