Bemelmans J. et. al. (editor) Elliptic and Parabolic Problems: Rolduc and Gaeta 2001

Подождите немного. Документ загружается.

410

[6] A. Hofler and N. Strecker, On the Coupled Diffusion of Dopants and Silicon Point

Defects, Technical Report

94/11,

Integrated Systems Laboratory, Swiss. Fed. Inst, of

Technology (ETH), Zurich, 1994.

[7] A. Hofler, Development and Application of a Model Hierarchy for Silicon Process Sim-

ulation, PhD thesis, Swiss. Fed. Inst, of Technology (ETH), Zurich, 1997. Also available

as Series in Microelectronics, vol. 69, Hartung-Gorre, Konstanz, 1997.

[8] Y. Saad,

ILUT:

A dual threshold incomplete ILU factorization, Numerical Linear Al-

gebra with Applications 1, pp. 387-402, 1994.

[9] A. Pomp, O. Schenk and W. Fichtner, An ILU Preconditioner Adapted to Diffusion

Processes in Semiconductors, Technical Report

99/11,

Integrated Systems Laboratory,

Swiss. Fed. Inst, of Technology (ETH), Zurich, 1999.

[10] S. Selberherr, Analysis and Simulation of Semiconductor Devices, Springer, Wien,

New York, 1984.

[11] P.A. Markowich, The Stationary Semiconductor Device Equations, Springer, Wien,

New York, 1986.

[12] P.A. Markowich, C.A. Ringhofer and C. Schmeiser, Semiconductor Equations,

Springer, Wien, New York, 1990.

[13] S.M.. Sze (ed.), Modern Semiconductor Device Physics, Wiley, New York, Chichester,

1998.

[14] X. Huang et al., Sub-50 nm P-Channel

FinFET,

IEEE Transactions on Electron

Devices 48, pp. 880-886, 2001.

Numerical Computation of Electromagnetic Guided

Waves in a General Perturbed Stratified Medium

Gomez Pedreira,

Dpto.

de Matematica Aplicada,

Univ. de Santiago de Compostela, Spain

Email: malola@usc.es

1 Introduction

The computation of electromagnetic waves propagating in perturbed stratified waveguides

is important both for the theoretical analysis and many applications in physics and engi-

neering, like integrated optics. An open stratified waveguide is a structure composed by

parallel layers of dielectric materials whose electromagnetic behaviour is characterized by

a physical property called the refraction index n.

The goal in the study is to compute the frequencies for which harmonic electromagnetic

waves of finite transverse energy can propagate without attenuation along the waveguide.

These waves are called guided waves (or guided modes) and they are particular solutions

of the Maxwell's equations in the form

E(x,x

,t) = E(x)e

-^\

H(x,x

,t) = B.(x)e^

px3

\ ^ '

where u > 0 is the angular frequency of the wave, ft > 0 is the wavenumber (or propagation

constant) of the mode, and E(x) = {Ei(x),

{x),

(x)) and ffl(x) = (Hi(x),

(x),

(x))

are vector functions of x = (11,2:2) describing the distribution of the electromagnetic field

in each cross section. Their main property is that their amplitude has finite transverse

energy

f (|E|

+ |M|

)da; < 00. (2)

In fact, it is this condition which determines if a mode is guided or not and physically

means that the energy of the mode remains practically confined in some bounded region

of the cross section.

We will assume that the guide is invariant, with respect to both the geometry and

its physical characteristics, under any translation in the propagation space direction x

Thus the refraction index will depend only on the two transverse coordinates (xi, x

The existence of the guided modes will depend on the refraction index distribution in

the cross section of the guide. If the materials composing the different layers are chosen

with a suitable refraction index, then the energy of the wave is vertically confined in the

411

412

material with the largest index. Nevertheless, to confine the light laterally —which is

essential for the design of these devices— it is necessary to modify the refraction index

in a bounded region, here denoted by K, in such a way that if the refraction index in

this region is properly chosen, the energy of the wave remains practically confined in a

neighborhood of

K..

The dimensions of

are so small in comparison with the stratified

medium that, from the mathematical point of view, the cross section of the guide is

considered as an unbounded domain. This is why this kind of guides are called open

waveguides.

We will assume that

n{x) = n(x

) Vi = (x

)^Kcl

where n £ L°°(R) is a function only of the vertical variable which represents the refraction

index of the so called reference medium, which is nothing but the stratified medium

without the perturbation K.

Mathematically, the problem to be solved can be reduced, under the assumption of

weak guidance (i.e. large wavenumber and weak variations between the maximun and the

minimun values of the refraction index, see Bonnet[l], Vassallo[2] ), to

(V)

Find w > 0,

> 0 and u € L

) (ti/ 0) such that

- A u +

u =

where u(x) = Ei(x),

{x),

Hi(x) or

(x).

For P > 0 given, (V) is a scalar eigenvalue problem set in the cross section of the waveg-

uide,

where w

plays the role of eigenvalue of the selfadjoint operator Ag = l/n

(—A + /3

)

of L

) and u, which denotes any of the transverse components of the electric or mag-

netic field, is the associated eigenvector.

We will be interested in the computation of guided modes associated to eigenvalues J

in the discrete spectrum of Ag, assuming they exist (see Bonnet et al. [3] or G6mez[4] for

existence conditions). These are isolated eigenvalues of finite multiplicity which satisfy

(cf. Gomez and Joly[5])

< u? < <j

(/3), with n+ = sup n(x), (3)

xeK

and <7

(/3) denotes the greatest lower bound of the essential spectrum of the operator Ag.

In the sequel, we shall denote

E = {(w, ffl6R

>0,^>0,

ln\

< aM }. (4)

Remark 1.1 Let us point out that the guided modes, even when they exist, do not exist

for just any values of

and

but for u and/3 linked by a certain relation

= f(/3) called

the dispersion relation of the mode.

From the numerical point of view, the main difficulties in solving the problem (V)

come from the unboundedness and the stratified nature of the domain, which made the

numerical methods proposed by other authors to be either non applicable , either "non

413

Figure 1: Sketch of an open stratified waveguide and the model cross section

exact" or very costly from the numerical point of view (cf. Joly and Poirier[6] or Mahe

[7])-

The main idea of the method proposed in this work for the computation of the guided

waves consists of a reformulation of the original problem (V) into a new one via the

introduction of a pseudo-differential operator K depending on (u,

j3)

in such a way that

solving the problem (V) is reduced to finding the values of (w,/3) for which —1 is an

eigenvalue of

K(LO,

/3).

The efficiency of the method relies in an efficient way to compute numerically this

operator. The practical evaluation of K(o;,

j3)

is obtained by a combination on ana-

lytical computations —which take into account the unboundedness of the propagation

medium— and numerical approximations involving: a domain truncation parameter R,

a series truncation parameter N and a space step h for the finite element approximation

of an auxiliary problem localized in a neighborhood of

K,.

The method combines three

different techniques: Fourier transform, Fourier series and mixed finite elements.

A particular case of this method was introduced and analyzed in Gomez and Joly[5, 8]

for the case of a medium composed of three homogeneous layers with a perturbation

K, of the refraction index involving only the central one. In these papers, the authors

have announced that such a method was generalizable to the case where the interior

layer is itself a stratified medium (made up of parallel layers), where the perturbation is

not necessarily embedded in one of them. The objective of this work is to present the

generalization of the method introduced in Gomez and Joly[5].

2 The new formulation of the problem

2.1 A first reduction to the boundary T

Let us going to consider a reference stratified medium composed of three different domains:

two exterior homogeneous layers —denoted by fl+ and f2~ — separated by an interior

medium

Q,,

= {(xi,X2) £R

,0<i2< L} made up of q parallel homogeneous layers. For

j = 1,..., q , we denote hj-\ and hj the x

coordinates corresponding to the lower and

upper boundaries respectively of the jth layer Oj (see Figure 1). Thus, the refraction

414

index of the reference medium is a piecewise constant function of x

given by

if x

€ £l~ = {(x

) e R

, x

< 0},

n(x

)

if x

€ Qj = {(2:1,2:2) G R

, < /ij_i < £

< /ij},

+ if i

G(l+ = {(ii,i

)eK

,i2>I},

(5)

where (5), rij denotes the constant value in the jth layer and we will assume, without loss

of generality, that n+, > n^,.

The first step of the method is to formulate a problem whose unknown

= (ip

<p~)

is the trace of u (the solution of (V)) on V = r

, with

r- = {(zi,0), Xi eR} = 8Q-,

r+ = {(xi,L), Xi 6 1} = 90+.

For (w, ft) £ E given, we introduce the operator S(ui,

defined as follows (for the sake of

simplicity, we omit to mention the appropriate functional framework and refer the reader

to Gomez and Joly[5, 8] for more details):

S(w,P)<p =

du(ip)

du{ip)

'du(tp)~

where

u(<p)

denotes the solution of the boundary value problem

)

-A

u(<p)

+ (P

- n

)

u(ip)

= 0

in Q = .

on T.

\r,

(6)

(7)

and n denotes the outer unit normal vector, with respect to the domain Q;, on boundaries

T+ orT". Moreover, [q]

{[q]r+,[q]r-),

[?]r+ = (9|n

)|r+ " (9|n+)lr+> Mr- = (9|

)|r- ~

(9|n.-)lr-

The idea is the following: take a function

defined on T and solve the boundary value

problem (V^), (which consists, in fact, of two decoupled problems, one outside the strip

Qi and another one inside). By construction, the function

u(tp)

is continuous. In order

that u{if) be a solution of the problem (V), it is enough to ensure the matching of the

normal derivatives on the lines x

= 0 and x

= L. This means that the jump of the

normal derivative of

u(tp)

across T, namely

S(u},/3)ip,

must be equal to 0. Then, problem

(V) is equivalent to

CPs)

For p > 0, find w > 0 with (w, p) 6 E

such that 0 is an eigenvalue of S(uj,

/3).

For numerical purposes, we use a decomposition of S{ui,P) into three operators

S(u,

= Si{u, p) + S

(w,

- S

(cj,

(8)

which corresponds to a decomposition of u as:

415

• u|

=uf, with uf = u*(tp) the unique solution of

f -Au± + (p-n^rf)ut = 0 in S£,

(Fe}

I «± = v on r,

• u|

. = Uj +

Up,

with Uj =

Ui(<p)

the unique solution of

-A

+ (fi

-n

)ui = 0 in n*,

l_ M; =

on T,

and u

= u

{ip) the unique solution, for certain values of a; and p, of

—A Up

+ (P

-n

= {n

-n

Ui in Q

Vpl

\ u

= 0 on T.

Then we define „, „, (dui

du~ \

(u,P)<p = —

~\ , —±-\

on lr+ on lr-

Si(u,P)<p

fdw 3M,| \

^ , > (dun

\ on |r+ on

Problems (

) and (T'i) are coercive (see Gomez and Joly[5]) and hence uniquely

solvable. In Gomez and Joly

[5],

it is also shown that problem (V

) is well-posed if and only

if u

g Gi(P), where G

(

8) is a finite (possibly empty) set of irregular frequencies (such

as those appearing in integral equations) which can be treated and computed separately.

Thus,

£ Gi(fi), S

, and thus S, are well defined.

2.2 Introduction of the operator K(w,/3)

We refer to Gomez and Joly[9, 5] for various properties of Si{ui,P), S

{u),P) and S

{ui,p).

In particular, it is explained that S(w,P) is not very easy to handle numerically, because

it has a continuous spectrum. That is why we introduce the operator

K(w,/3) = (S

-S

, (10)

which is well defined as a linear operator in L

(D since Si

—

is an isomorphism from

H^r) into L

(r). It can also be proven that K(w,

is a compact operator whose spectrum

is purely discrete and consists of a countable infinite set of real eigenvalues admitting 0 as

unique accumulation point (the proof is based on the fact that S

is a compact selfadjoint

operator in L

(D, see Gomez and Joly[9]).

The relationship between

S{CJ,

and K(u,

is simply (I denotes the identity operator)

S(co,p) = (S

-S

){I + K(u

,p)}, (11)

so it is obvious that the problem (Ps) is equivalent to the following

For a given

/3,

find w > 0 (w

0 Gj(/3)) with (w,0) e E,

)

' such that -1 is an eigenvalue of

K(LJ,P).

416

3 Computation

the operators

, Si

and

The practical computation of the operator

K(OJ,/3)

involves the computation

the op-

erators 5

and

The main difference with the method introduced

Gomez and

Joly[5] lies in the computation of the operators

and

Let us briefly explain

this

section how these computations are done (see Gomez and Joly [10] for details).

Exploiting the fact that both (V

) and (Pi) are invariant under translation

the

direction, these problems can be solved by partial Fourier transform in

variable. As

consequence, the operators

and Si, and then operator (Si

—

)~\ can be analytically

computed.

Let us denote by

the dual variable of x\. By solving the problem (P

) in the Fourier

domain and applying the definition of the operator

given in (9),

is easy to check that

[Sefw^te] (*)

(u,

ftk)

where

J & =

(fc

+ ^

-n+^)i/2

and

When solving the problem (Vi) we must take into account that the domain Q;

is a

stratified medium composed by q parallel layers. The equation (in the Fourier domain)

-xi-

(fc

+ p

n}ui

)

= 0, j =

1, • • •

, q

satisfiyed

each layer must

those

the adjoining layers

by the

continuity of the solution

and of its normal derivative.

The computations lead

[SiK/3")?] (*) = M

{

(uj,P; k)

ftk).

where the entries

the matrix Mi(uj,P;k) associated

the operator

Si in

the Fourier

domain are given by

{LJ, P;

k) = 6 (e-^

+ T

) + T

+ T

)/D

(u,P;k)= -2^e~^

(uj,P;k)

-2^e~^

(w,

£, (T

^(T

))/D

with

D =

-T

417

and [T(

](,

n=lj2

denotes the elements of the matrix T which results from the product of

matrices T

_! T

• • •

7\ given by

T-

-. (

hj (fi+l-fj) (1

JL.)

-hj fe+i+fj) (i _ Ji_\

Matrices Tj, j = 1,... ,q

—

1 are called transmission matrices since they derive from the

application of the boundary conditions to be met at an interface between two layers.

Remark 3.1 The fact that D ^ 0 for any k is not immediate but is a consequence of the

fact that problem (Vi) is coercive.

Contrary to the operators S

and Si, operator S

cannot be analytically computed.

The definition of this operator involves the resolution of the problem (V

) which is a

boundary value problem posed in the infinite strip 0*. Taking into account that the

right hand side has compact support included in /C and the fact that the solution u

satisfies and homogeneous Dirichlet condition, this problem can be numerically handled

by using a localized finite element

method.

This method will allow us to reduce the

actual computations to a rectangle Of, surrounding the perturbation with the help of a

transparent boundary condition. The method implies introducing two artificial vertical

boundaries in the domain Q

namely S

= {a

} x (0,L) and E~ = {a

} x (0,L), to

split

Cli

into a bounded domain 0(, containing the perturbation and another one Qj \ 0(,

where the refraction index is a piecewise constant function of x

(see Fig. 2). The boundary

conditions to be defined on E* involve the introduction of a Neumann to Dirichlet operator

T(ui,/3) depending on w and

given by

[T(w, flip]{x

) = 53 CV fl w

)

(13)

*=i

where

• £fc(w, P) =

(/3

+ Ak(w))

, which is strictly positive as a consequence of the inequal-

ity (3).

• Afc(ai) denotes the eigenvalues of the one-dimensional eigenvalue problem

-w"(x

) -n

)u}

w(x

) = Aw(i

), w(x

)

w(0) = w(L) = 0.

[

• {wk(u},x

)} denotes the set of associated eigenfuntions, which forms an hilbertian

basis of L

(0,L).

•

tpk

denotes the expansion coefficients of

in the basis

{w^w,:^)}-

By using the operator

T(CJ,

ft), we can write an equivalent (in a certain sense) formula-

tion of the problem (V

) which couples a variational formulation in the bounded domain

Qj with a Fourier expansion outside. More precisely, the computation of u

has to be

done in two steps:

418

One computes u

inside

by solving numerically the boundary value problem

_!_ fl

— m

/.,

\ ti — In

— m

\ ,.,

3£\

re)

—Aup + (^

—

) u

= (n

—

) ai

«*

in fi|,,

u„ = o on r

= r n dn

-T(u,,P)^= u

on E±,

where u denotes the outgoing unit normal vector, with respect to the domain Q

on boundaries E

or S~.

• Knowing u

inside Q

, it is analytically computed in the exterior domain Q

\f2(, via

the formulas

) = ^2u^w

{u},x

)e-

ik{xi

a+)

if xi > o+,

(15)

) = ^2,ulw

{u,x

)e~

ik{a

x,)

if xi < a~,

where u£ and u^ denote the expansion coefficients in the basis {w^} of the trace on

E+ and £~ of the solution u

of (Vf).

Once u

has been computed, then it suffices to apply the definition (9) to compute S

4 Numerical approximation of the problem (VK)

4.1 The truncation of

The first difficulty one has to face for solving problem (VK) lies in the fact that the

unknown function tp is defined on a one dimensional but unbounded domain, namely

the boundary T. To avoid this problem, Gomez and Joly[8] propose an approach based

on a truncation method that can be applied for a multi-layered structure like the one

considered in this work following an analogous procedure.

The idea is working with functions defined on the bounded domain

= T n {(x

) I a~ - R < xx < a+ + R},

where R > 0 is an approximation parameter which goes to +oo.

Let us consider the orthogonal projector

IIR

on L

) defined by

: L

(r) —• L

(r)

—•

<p,

(16)

where x

denotes the characteristic function of T

. The idea is to write an equation for

TlRtp,

the "restriction" of

to T

. Such an equation does not exist but we can write an

approximate equation, up to exponentially small errors, whose unknown tp

will be an

approximation of

TlRtp.

419

~" I Q„

l!I

S+

series expansion

Figure 2: Numerical approximation

Error estimates that can be proved thanks to regularity and decay properties of the

solution Up of the interior problem (V

) (cf. Gomez and Joly[8]), allow us to approximate

the problem {PK) by

KRK

^eL

(r„)

(

(17)

with

(Si

(is)

This method has been proven to be exponentially accurate (cf. Gomez and Joly[8]). More

precisely, it can be proven that the nonzero eigenvalues of K

converge exponentially to

the nonzero eigenvalues of K as R goes to infinity. In particular, the eigenvalue —1 will

be approximated with an exponential accuracy.

4.2 The series truncation

In order to solve numerically the problem {Vp) we must truncate both the series in (13)

and in (15) at (large) rank JV. Then, the problem (Vp) is approximated by

-Att£ + (/?

- n

)w£ = (n

- n

)uj

in Q

<= 0

- u"

on T

on E±

(Vf)

{

where T

denotes the operator obtained by the truncation at rank N of the series given

by (13).

Once time we have computed u

in Q

, the solution u

e fij \ Oj, is computed via a

Fourier series development

*=i (19)

) = '^2

(k{a

Xl)

ifa:i<o"-

Taking into account the series truncation, the operator S

is approximated by the operator