CONTENTS
transversally limited with various boundary conditions. The mode orthogo-
nality is analytically proven when the anisotropy is only dielectric (or mag-
netic) and confined in a plane normal to the waveguide axis.
Successively, a first extension of the standard SEM is introduced. In
Chapter 4 the problem of finding guided modes in open planar structures
is tackled and solved using suitable expansion functions in the unbounded
regions in order to model the damped fields.
Finally, in Chapter 5 the standard SEM is applied to find modes in cylin-
drical layered waveguides comprising isotropic media, after writing two wave
equations for the longitudinal fields.
Every chapter ends presenting a study of the convergence properties of
the SEM; in particular the exponential rate of convergence is justified in all
practical cases and normalized error plots, which are useful in evaluating and
predicting the error committed on the approximated solutions, are obtained.
Numerical results such as dispersion curves and field profiles are also
reported as some example of application.
Two Appendices complete the work, providing the main formulae rela-
tive to Legendre and Laguerre polynomials and the matrices that allow to
implement the SEM automatically.
2
Chapter 1
Introduction
1.1 Applications of layered media
Many optoelectronic devices and components such as directional couplers,
modulators, semiconductor lasers and optical interconnections employ biaxial
media, both natural and artificially generated by electro-optic effect.
Since the device can often be modeled as a planar layered waveguide,
the design and optimization of these structures require accurate and efficient
methods to compute the propagation characteristics of the guided modes.
Also multilayered elastic media support the propagation of various types
of waves and the determination of the field functions and propagation con-
stants is a classical problem. Finally, optoelectronics uses devices in which
acoustic and electromagnetic waves interact.
The eigenmodes can be utilized to express the relevant Green’s function,
which is necessary to formulate scattering and junction problems by the
integral equation technique (see for example the works [1], [2]).
Electromagnetic planar stratified media are generally analyzed by char-
acterizing each layer by a transmission matrix [3]-[6] or a scattering matrix.
Even the elastic media are typically analyzed along the direction normal
to the interfaces: the T-matrix [7], [8], the propagator matrix [9], coupled
transmission lines and lumped circuits [10] are commonly employed in de-
scribing the layers.
3
1 – Introduction
Then, each elementary block is connected to the remaining ones by stan-
dard cascading formulas. The propagation constants β of the structure can
be found by imposing the well known transverse resonance condition, i.e.
searching for field configurations that can exist in absence of sources. This
procedure leads to a transcendental equation in the variable β.
Another typical method of studying layered media lays on the fact that
plane waves can propagate in an unbounded medium and thus they can be
seen and used as a basis to represent the field in a longitudinally limited
medium.
In practice the field is expanded by superposing four (six) partial plane
waves in the electromagnetic (elastic) case. The amplitude of each wave is
unknown and can be suitably determined in order to achieve the continuity of
some field components. These conditions lead to a homogeneous system; non
trivial solutions exist only for the particular values of β that reduce to zero
the determinant of the coefficient matrix [12]. Again the resulting equation
is transcendental.
In either case analytic function theory may be applied to locate the roots
in the Gauss’ plane, however automated search still remain a rather difficult
task, in particular when losses are present, and some solutions may be missed.
An interesting technique that assures the capture of all bound modes
for lossless and open waveguides is presented in [14], but it is restricted to
isotropic layers.
To overcome the above mentioned disadvantages, a numerical method to
determine mode functions and propagation constants of a layered structure
without solving any transcendental equation is proposed and discussed. The
method is based on the representation of the fields in terms of basis functions
and the consequent reduction of the differential eigenvalue problem to an
algebraic one.
A similar technique was employed in [15] to study optical gratings com-
prising isotropic dielectric layers.
In principle, an entire domain basis might be chosen, (see, for exam-
ple, [16], [17], where Laguerre polynomials were used to investigate acoustic
surface waves in semi-infinite layered media), but this would cause a slow
(power-law) convergence, owing to the fact the field component are infinitely
4
1 – Introduction
differentiable everywhere in space but at an interface.
Hence, separate bases in each layer are required to ensure a correct mod-
eling of the discontinuities. An exponential convergence ensues adopting
Legendre polynomials as expansion functions.
Both in acoustic and in electromagnetic structures the propagation phe-
nomena are described by second order differential equations. Therefore in
the following sections, regardless to the physical nature of the problem, the
method just outlined will be examined in details and applied to an abstract
differential problem.
1.2 Abstract formulation for general linear
fields
A somewhat general differential equation (at least for the purposes of the
present work) can be written formally as
L
(
x,
∂
∂x
; ζ(x)
)
ϕ(x) = λM
(
x,
∂
∂x
; ζ(x)
)
ϕ(x) (1.1)
where the independent variable x ranges in the closed interval [a,b]. Succes-
sively, it will be clear that many actual physical problem can be formulated
in the form (1.1).
L and M represent two linear differential operators that may depend on
∂/∂x and x both directly and through ζ, which in turn stands for a set of
parameters pertaining to the physical problem (e.g. the permittivity and the
permeability in electromagnetics, the stiffness and the density in acoustics
and so on). The scalar quantity λ is called the eigenvalue, while ϕ(x) is the
eigenfunction.
Though it is not necessary for the successive considerations, however we
will suppose at most the presence of second order derivatives in L and M,
since a wide class of physical problems can be modeled by second order partial
differential equations. As it will appear further on, this choice will not affect
the description of the spectral element method.
As is well known, the general integral of Eq. (1.1) can be written as a
linear superposition of two independent solutions (on assuming at most the
5
1 – Introduction
presence of second order derivatives) by means of two constants that are
arbitrary at all. In addition, λ has indeed no restriction and may take any
value in the complex plane. On the other hand, from a physical point of view
only those solutions satisfying certain conditions at the ends of the interval
[a,b] can be accepted and properly interpreted.
So we are concerned more precisely with the solution of Eq. (1.1) supplied
with two boundary conditions we can formally state as
B
(
x,
∂
∂x
; ζ(x),λ
)
ϕ(x)
∣
∣
∣
∣
∣
x=a,b
= 0. (1.2)
These two constraints force λ to assume only an infinite but discrete number
of well determined values. The problem can be reformulated like this: finding
all the possible pairs {λ, ϕ(x)} that solve Eq. (1.1) and satisfy Eqs. (1.2).
Something more has to be said on the form and nature of the eigenfunction
ϕ. To be more specific, let H be the Hilbert space of the complex valued
functions of real variable x, ranging in the interval [a,b],
f : [a,b] −→ C (1.3)
with inner product
〈f,g〉 =
∫ b
a
f(x)g∗(x)dx. (1.4)
Then, we assume ϕ(x) to be a column vector with Q function elements be-
longing to H, i.e.
ϕ(x) =
ϕ1(x)
.
.
.
ϕQ(x)
, (1.5)
so that ϕ(x) can be considered as an element of the Hilbert space HQ =
H×H× . . .×H, with inner product
〈ϕ,θ〉 =
Q∑
q=1
〈ϕq,θq〉 =
Q∑
q=1
∫ b
a
ϕq(x)θ∗q(x)dx. (1.6)
Though in principle, if ζ(x) is continuous, ϕ(x) is C∞ in [a,b], we suppose
that it is not differentiable in a finite number of points, the reason of which
will be made clear further on.
6
1 – Introduction
L, M and B are operators that, acting on an element of HQ, transform
it into an element of the same space. Therefore, they can be represented by
symbolic matrices with size Q×Q, e.g.
L =
L11 · · · L1Q
.
.
.
. . .
.
.
.
LQ1 · · · LQQ
. (1.7)
At this point it is worth mentioning that a bit more general kind of dif-
ferential equations, called polynomial eigenvalue problems, where λ appears
raised to more than one integer power, can be reduced to Eq. (1.1) and then
treated in the same manner.
As a practical example consider
A0ψ(x) + λA1ψ(x) + λ2A2ψ(x) = 0. (1.8)
After introducing an auxiliary unknown function, say
τ(x) = λψ(x), (1.9)
(note that the presence of the eigenvalue in defining τ(x) is a basic feature
of this reduction process), Eq. (1.8) can be restated as follows
A0 A1
0 I
ψ
τ
= λ
0 −A2
I 0
ψ
τ
, (1.10)
with I representing the identity operator; now it is very easy to identify
L, M and ϕ. Other choices are also allowed; in fact, L and M might be
constructed in a different way and moreover the auxiliary unknown τ might
be chosen as λL2ψ.
To proceed further, we will limit our study to the problems with ζ(x)
being a piecewise constant function of x; in particular, we will assume that ζ
takes N distinct values in as many subintervals of [a,b] = [x1,xN+1], namely:
ζ(x) = ζ i, for x ∈ ]xi,xi+1[, i = 1,2, . . . ,N. (1.11)
The values ζ(x) takes on in the points xi are not relevant at all. In the i-th
subinterval Eq. (1.1) becomes
L
(
x,
∂
∂x
; ζ i
)
ϕi(x) = λM
(
x,
∂
∂x
; ζ i
)
ϕi(x) (1.12)
7
1 – Introduction
or more succinctly,
Liϕi = λMiϕi (1.13)
where the operators Li and Mi are defined only on ]xi,xi+1[.
The N differential equations represented by (1.13) could be solved sep-
arately and consequently N solutions, each holding in the i-th subinterval,
would be found. Of course, it is also clear that this new differential problem
is not perfectly equivalent to the original stated by Eqs. (1.1) and (1.2). In
fact, if no more condition is imposed, once again λ may assume arbitrary
values.
However, as Eqs. (1.13) and (1.2) actually model a physical phenomenon,
more informations about the behavior of the whole solution ϕ(x) are surely
available, and in general they just specify the way the ϕi must be connected.
These further conditions turn into 2(N−1)Q constraints at the points xi,
since we have supposed the existence of second derivatives in the operators.
As it is obvious that in the points where ζ(x) is discontinuous ϕ may exhibit a
singular behavior, this justifies the assumption made above on the differential
properties of ϕ(x).
Regarding the kind of conditions imposed on the eigenfunction in presence
of discontinuities of ζ(x), the continuity of ϕ(x) is generally required:
ϕi(xi+1) = ϕi+1(xi+1), i = 1,2, . . . ,N − 1 (1.14)
and these Q(N − 1) equations imply that ϕ is at least C0 in [x1,xN+1]. The
remaining needed conditions can be formulated as follows:
S
(
x,
∂
∂x
; ζ i
)
ϕ(x)
∣
∣
∣
∣
∣xi+1
= S
(
x,
∂
∂x
; ζ i+1
)
ϕ(x)
∣
∣
∣
∣
∣xi+1
(1.15)
Sometimes, for a certain class of problems, these last equations also assure
that ϕ(x) is C1 in [x1,xN+1], but in general this property does not hold.
The problem is again well posed: we have to find all the possible values
of λ and the corresponding functions ϕi(x) that solve Eq. (1.1), subjected to
conditions (1.2), (1.14) and (1.15). The eigenvalue λ can assume an infinite
but discrete number of values and the related eigenfunctions ϕ(x) constitute
a basis of the space HQ.
8
1 – Introduction
In certain cases, the orthogonality of the eigenfunctions, in the sense of
the inner product (1.6), can be also proved; if they were not orthogonal, the
adjoint problem should be solved, as the orthogonality to the adjoint basis
is guaranteed.
1.3 The Spectral Element Method (SEM)
Equation (1.1) supplemented with the boundary conditions (1.2) could be
solved directly by the moment method, choosing the sets of basis and test
functions both C∞ on [a,b]. However, it was outlined in the previous Section
that ϕ(x) can be differentiated only once. Hence, representing the solution
by means of functions that are infinitely differentiable (typically sines and
cosines or complex exponentials, that in addition are periodic) does not con-
stitute the best choice. As a consequence, a Gibbs phenomenon may be
expected as well as a slow convergence of the approximate solution to the
exact one.
In order to overcome these two undesired effects, a valid and effective
alternative consists in expanding the unknown separately in each subdomain,
using different sets of functions. In other words, it is more convenient to face
the solution of Eqs. (1.13). The various conditions (1.2), (1.14) and (1.15)
are enforced directly, so that they are always fulfilled by the approximate
solution, though the sets of expansion functions (theoretically infinite) are
truncated to a certain point for computational purposes.
As the number of subintervals is fixed (very often by the nature of the
physical problem under study) and the accuracy of the approximate solution
is increased inserting a larger and larger number of expansion functions in
each subdomain, the technique we are going to describe belongs to the class
of spectral element methods. Furthermore, if the basis functions do not satisfy
individually the various constraints in the nodal points xi, the technique is
also called a tau method and is due to Lanczos [18].
In principle, any set of functions infinitely differentiable on a bounded
interval might be used to represent the unknowns ϕi. A power series seems
very simple to handle, but the basis {1,x,x2, . . .} is not orthogonal and this
fact would give rise to numerical instability. We will adopt bases of Legendre
9
1 – Introduction
polynomials, which are orthogonal (see Appendix A for details) in the interval
[−1,1] and can be quickly calculated in a numerically stable manner by means
of the recurrence formula (A.2).
Since the Legendre polynomials are orthogonal in [−1,1], but we have to
expand functions defined on a generic interval of the x-axis, N changes of
independent variable are needed in order to map each subdomain [xi,xi+1]
onto the orthogonality interval:
ξi = αi (x− xi) (1.16)
with
αi =
2
xi+1 − xi
, (1.17)
xi =
xi+1 + x1
2
. (1.18)
Our aim is now to reduce the differential problem described in Section
1.2 to an algebraic generalized eigenvalue one, whose solution is a classical
and well known question of numerical analysis. To do so, we define N bases
as
{Pm(ξi)}Mim=0, i = 1,2, . . . ,N, (1.19)
and then expand the unknowns:
ϕi(ξi) =
Mi∑
m=0
Pm(ξi) pim,1
.
.
.
Pm(ξi) pim,Q
. (1.20)
Substituting these expressions in Eqs. (1.13) yields
Mi∑
m=0
Li
Pm(ξi) pim,1
.
.
.
Pm(ξi) pim,Q
= λ
Mi∑
m=0
Mi
Pm(ξi) pim,1
.
.
.
Pm(ξi) pim,Q
. (1.21)
Since now the independent variable is ξi, obviously some modifications occur
in the elements constituting the operators, but we postpone the details to
the next section.
10
1 – Introduction
These N equations are projected, according to the inner product defined
by Eq. (1.6) and restricted to the interval [−1,1], not onto the bases (1.19)
as expected, but onto the following ones:
{Pm(ξi)}Mi−2m=0 , i = 1,2, . . . ,N, (1.22)
and the explanation of this is as follows.
Let us suppose for the moment to perform a projection on the bases
(1.19). A generalized eigenvalue problem would follow immediately, since the
two matrices appearing in the two sides of the equation would be intrinsically
square, but it becomes also evident that the various conditions at the nodal
points have not been employed.
At this point one can observe that these constraint equations are stronger,
as they guarantee the proper physical behavior of the searched solution.
Thus, it is more convenient to force them in place of the last equations
following from the differential system.
However, there is perhaps a more important feature regarding the projec-
tion process. The last two equations following from each of (1.21) are linearly
dependent and this is strictly related to the fact that the second derivative
of a polynomial with degree Mi has degree Mi − 2. It will be seen that in
certain cases the above mentioned equations reduce to
c
p
iMi−1,q
piMi,q
= λ
p
iMi−1,q
piMi,q
. (1.23)
They are uncoupled from the others and, since in general the eigenvalue
differs from the scalar c, they lead to the trivial solution piMi−1,q = piMi,q = 0.
This observation further justifies the use of the constraint equations.
In practice, a projection onto only the first Mi− 1 elements of each basis
set (1.19) is performed, thus yielding
Q∑
q=1
Mi∑
m=0
⟨
Pn(ξi),LirqPm(ξi)
⟩
pim,q = λ
Q∑
q=1
Mi∑
m=0
⟨
Pn(ξi),MirqPm(ξi)
⟩
pim,q, (1.24)
where n = 0, . . . ,Mi − 2 and r = 1, . . . ,Q. These represent Q∑Ni=1(Mi − 1)
equations in the Q∑N
i=1(Mi + 1) unknown expansion coefficients pim,q.
The constraints provide the missing relationships to complete the alge-
braic system, in particular:
11
1 – Introduction
• Q equations from boundary condition in x = x1 (or ξ1 = −1)
M1∑
m=0
B
(
x,
∂
∂x
; ζ1,λ
)
Pm(ξ1) p1m,1
.
.
.
Pm(ξ1) p1m,Q
∣
∣
∣
∣
∣
∣
∣
∣
∣
ξ1=−1
= 0, (1.25)
• Q equations from boundary condition in x = xN+1 (or ξN = 1)
MN∑
m=0
B
(
x,
∂
∂x
; ζN ,λ
)
Pm(ξN) pNm,1
.
.
.
Pm(ξN) pim,Q
∣
∣
∣
∣
∣
∣
∣
∣
∣
ξN=1
= 0, (1.26)
• Q(N − 1) equations from continuity conditions in x = xi+1, with i =
1, . . . ,N − 1
Mi∑
m=0
Pm(1) pim,1
.
.
.
Pm(1) pim,Q
=
Mi+1∑
m=0
Pm(−1) pi+1m,1
.
.
.
Pm(−1) pi+1m,Q
(1.27)
• Q(N − 1) equations from further conditions in x = xi+1
Mi∑
m=0
S
(
x,
∂
∂x
; ζ i,λ
)
Pm(ξi) pim,1
.
.
.
Pm(ξi) pim,Q
∣
∣
∣
∣
∣
∣
∣
∣
∣
ξi=1
=
Mi+1∑
m=0
S
(
x,
∂
∂x
; ζ i+1,λ
)
Pm(ξi+1) pi+1m,1
.
.
.
Pm(ξi+1) pi+1m,Q
∣
∣
∣
∣
∣
∣
∣
∣
∣
ξi+1=−1
(1.28)
All these equations together with Eq. (1.24) can be written in a compact
form:
[L][Φ] = λ[M ][Φ]. (1.29)
The order of the equations is rather arbitrary and it can be exploited to give
the matrices a simple structure more suitable for computational purposes.
On the other hand, the columns of the matrices [L] and [M ] are strictly
12
1 – Introduction
related to the way of ordering the unknown coefficients collected in [Φ]. A
suitable choice seems the following one
[Φ] =
p1
1
.
.
.
pN
1
.
.
.
p1Q
.
.
.
pNQ
, (1.30)
where
pi
q
=
pi0,q
.
.
.
piMi,q
. (1.31)
In certain cases, if the boundary conditions do not involve the eigenvalue
λ, Eqs. (1.25-1.28) can be solved with respect to 2QN expansion coefficients
that are to be substituted in (1.24) to obtain a reduced homogeneous system,
say:
[L′][Φ′] = λ[M ′][Φ′]. (1.32)
This elimination process could seem a further complication, but in practice
it turns in time saving as the new eigenvalue problem has smaller size.
1.4 Construction of the algebraic system ma-
trices
In this section we want to show how the system matrices [L] and [M ] can be
constructed automatically by superposing few elementary blocks, provided
both the equations and the unknown coefficients are ordered in a suitable
manner. Their structure depends upon the particular basis functions chosen
(i.e. the Legendre polynomials), whereas the simple rules they obey are
related to the form of the differential operator.
Recall that the linear operators L and M are represented by a Q × Q
symbolic matrix as in Eq. (1.7). Moreover, their generic element, say Lrq,
13
1 – Introduction
takes N different values in each subdomain [xi,xi+1]. All these possible values
are inserted in a new formal operator as follows:
L˜ =
L111
. . .
LN11
· · ·
L11Q
. . .
LN1Q
.
.
.
. . .
.
.
.
L1Q1
. . .
LNQ1
· · ·
L1QQ
. . .
LNQQ
︸ ︷︷ ︸
Q block columns
. (1.33)
The form of [L] and [M ] will be now examined in more depth with the help
of the pictorial representation reported in Fig. 1.1. The following comments
apply indifferently to either the operators, but for simplicity we will refer to
[L].
The main Q columns correspond to the Q unknown functions, contained
in the vector ϕ. The first Q main rows originate directly from projecting the
field equation onto the reduced Legendre bases. Each block has in turn a
diagonal structure; the i-th element corresponds to the equations holding in
the i-th subinterval and has Mi − 1 rows and Mi + 1 columns. Therefore,
this part of the operator [L] has globally size Q∑i(Mi− 1)×Q∑i(Mi +1).
Coming down the rows of [L], the contribution of the continuity conditions
is encountered. A first stack of Q block rows comes from the continuity
of the unknowns, as stated by Eqs. (1.27); this part of [L] must result
block diagonal, since the unknowns ϕq are independent of each other. The
submatrices are band diagonal, as each of Eqs. (1.27) involves always a pair
of unknowns ϕiq in two adjacent subdomains. Similar considerations apply to
the second stack of Q block rows, that originate from Eqs. (1.28). Due to the
arbitrary and sometimes complicated nature of the operators S i, that may
couple the unknowns ϕq, this part of [L] is not block diagonal. However, the
submatrices are quite analogous to those following from the other continuity
conditions and this fact has been outlined in Fig. 1.1. This part of [L]
contains 2Q(N − 1) equations globally.
14
1 – Introduction
Block diagonal matrix
of field equations
0
01 2
N-1
N
Block band matrix of
continuity conditions
1 2
2 3
N-2 N-1
N-1 N0
0
Block submatrix of
boundary conditions
1 N0 0
11
Q1
1Q
QQ
11
Q1
1Q
QQ
11
QQ0
0
Figure 1.1. Schematic representation of the matrices [L] and [M ], showing
their block structure.
Finally, in the lower part of [L] the contribution of the boundary con-
straints takes place. Equations (1.25-1.26) may couple the unknowns and
they originate 2Q rows (Q at each external end). Of course, the form of the
submatrices depends on the nature of B
1
and B
N
; if they express either a
periodicity or pseudo-periodicity condition, they involve φ
1
q
and φ
N
q
and the
submatrix appears as in Fig. 1.1. In the remaining cases, they do not couple
the unknowns in the first and the last subdomain.
15
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.
