Technische Universität Berlin
Institut für Mathematik
On adaptive finite element methods for the simulation
of two-dimensional photonic crystals
Marine Froidevaux
Preprint 06-2018
Preprint-Reihe des Instituts für Mathematik
Technische Universität Berlin
http://www.math.tu-berlin.de/preprints
Preprint 06-2018 July 2018
On adaptive finite element methods for the simulation of two-dimensional
photonic crystals∗
MARINE FROIDEVAUX‡
Abstract. The first part of this paper is devoted to the modeling of wave propagation
inside a perfect two-dimensional photonic crystal and the spectral analysis of the resulting
eigenvalue problem. In the second part of the paper, we focus on a special case, where
the eigenvalue problem is linear and Hermitian. We introduce a residual-based estimator
approximating the algebraic error, i.e., the error in the computed eigenvector induced
by iterative methods. The estimator for the algebraic error approximates the error in
the same norm as the one used in already existing estimators for the discretization error,
therefore enabling error balancing between both types of errors, and thus improving the
efficiency of the adaptive finite element procedure.
Key words. Photonic crystals, error estimators, adaptive finite element methods
AMS subject classifications. 35P15,78M10,65F15
1. Introduction
Photonic crystals are composite materials, made of periodically arranged components,
whose periodic structure affects the propagation of electromagnetic waves. These waves
may propagate through or be reflected by the crystal depending on their frequencies ωand
on their wave-vector k, but also on parameters such as the electrical permittivity of the
constituent materials, or the geometry of the spatial arrangement of the crystal.
Photonic crystals creating so-called band gaps are of special interest for applications
in engineering. A band gap denotes a range of frequencies at which no wave can propa-
gate in any direction through the crystal. The design of photonic crystals for industrial
applications requires the optimization of many geometrical and material parameters, in
order to obtain large bad-gaps around certain frequencies specified by the application. For
this reason, we are interested in developing accurate mathematical models and numerical
methods which can enable the time-efficient simulation of these physical systems.
The mathematical modeling of a photonic crystal is based on the Maxwell equations. In
the case of a perfect infinite crystal, the Maxwell equation lead to a sequence of parameter-
dependant PDE eigenvalue problems (EVPs), which can be discretized via a conforming
finite element method. One can show that eigenpairs of the discretized problem converge to
the eigenpairs of the PDE EVP as one increases the number of basis functions in the finite
element space, see e.g. [Hac92]. However, having time-efficient computations in mind, the
number of basis functions in the finite element space, as well as the spatial localization,
should be chosen in a way that minimizes the size of the discretized system, while warranty-
ing a certain accuracy. This kind of trade-off problems is typically addressed via adaptive
‡Institut für Mathematik, MA 4-5, Technische Universität Berlin, Straße des 17. Juni
136, 10623 Berlin, Germany
∗Research supported by the Einstein Center ECMath via project OT10 Model Reduction for Nonlinear
Parameter-Dependent Eigenvalue Problems in Photonic Crystals.
1
2 M. FROIDEVAUX
finite element methods, and error estimators are the central tool enabling this adaptiv-
ity. In [GG12], a reliable and efficient error estimator was developed in order to estimate
the error produced by the finite element method. However, for large-scaled problems it
is computationally too expensive to solve the discrete EVPs exactly via direct methods.
Instead, iterative methods are used, such as e.g. Lanczos’ method, in order to solve the
large-scale discrete problems up to a prescribed tolerance. In this setting, the error arising
from the solution of discrete EVPs via iterative methods, also called the algebraic error,
cannot be ignored as a contribution to the total numerical error. Moreover, as part of the
adaptive process, it may be sufficient to solve the coarsely discrete EVPs with an accuracy
of the magnitude of the discretization error. For these reasons, error estimators are not
only needed to approximate the discretization error, but also to approximate the algebraic
error. Moreover, both types of error estimators should be part of a unified framework,
where all errors are measured in the same norm, so that meaningful comparisons can be
made as part of error balancing procedures. For the Laplace eigenvalue problem, error
balancing between the discretization and the algebraic errors was addressed in [CDM+17]
and a fully-adpative algorithm called AFEMLA algorithm was introduced in [MM11], and
extensively discussed in [Mię10].
The mathematical description and simulation of photonic crystals has already been the
subject of many research papers. Discretization via Nédélec finite elements for three-
dimensional photonic crystals has been considered in [BCG06], while the same prob-
lem was solved in [HLM16] via a finite difference scheme followed by a Newton-type
numerical method. For two-dimensional photonic crystals, a spectral analysis was per-
formed in [GG12] for the special case where the electric permittivity is real and frequency-
independent. The extensive spectral analysis of the non-Hermitian non-linear case can be
found in [Eng10] and [ELT17]. Moreover, a linearization technique for the Drude-Lorentz
model of the electrical permittivity was introduced in [EKE12].
In the three first sections of this paper, we review the modeling of a perfect photonic
crystals based on the Maxwell equations which leads to the operator NLEVP. In the fourth
section, we consider the case where the electric permittivity is a real-valued and frequency-
independent function. We introduce the weak formulation of the resulting Hermitian linear
EVP and analyse its properties. In particular we show that its spectrum consists in at
most countably many real and positive eigenvalues with finite geometric multiplicities.
In the fifth section, we develop a residual-based error estimator which approximate the
algebraic error on the computed eigenvector in the norm defined by the sesquilinear forms
of the weak formulation. Finally, the last section of this paper is devoted to a numerical
experiment, where we simulate a two-dimensional photonic crystal and test the reliability,
as well as the efficiency, of the algebraic error estimator.
2. Preliminaries on the Maxwell equations
2.1. General Form. The Maxwell equations in their general form and in SI units are
stated as follows:
−∂tD(x, t) + ∇×H(x, t) = J(x, t),(Ampère’s circuital law)
∇·D(x, t) = ρ(x, t),(Gauss’ law)
∂tB(x, t) + ∇×E(x, t)=0,(Faraday’s law of induction)
∇·B(x, t)=0,(Gauss’ law for magnetism)
where E, H, D, B and Jare three-dimensional vector fields and ρis a real scalar field,
see e.g. [Jac99]. These fundamental equations come along with the following empirical
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 3
constitutive equations:
D(x, t) = 0E(x, t) + 0χE(x, t)∗E(x, t),
B(x, t) = µ0H(x, t) + µ0χM(x, t)∗H(x, t),
where 0, µ0are strictly positive constants, χE, χMare causal transfer functions, and ∗
denotes a convolution in time, i.e.
χE(x, t)∗E(x, t) =
∞
Z
−∞
χE(x, t −τ)E(x, τ)dτcausality
=
t
Z
−∞
χE(x, t −τ)E(x, τ)dτ.
Remark 2.1.In the remainder, we omit writing the explicit dependence of vector field on
space, time and frequency when the dependence is clear from the context.
In the remainder, we will assume that no free current nor free charges are present in
the system, i.e. J≡0and ρ≡0, and that the materials are non magnetizable, i.e.
χM≡0. Under these assumptions, the Maxwell equations simplify to the following system
of equations:
(2.1) −∂t(0E+0χE∗E) + ∇×H= 0, ∂t(µ0H) + ∇×E= 0,
∇·(0E+0χE∗E) = 0,∇·(µ0H)=0.
2.2. The eigenvalue problem. We recall briefly the definition of the Fourier transform.
Given a function f(t)∈L1(R), i.e., f(t)is absolutely integrable on R, its Fourier transform
is defined by
F{f}(ω):=ZR
f(t)e−iωt dt.
Note that when the variable trepresents time with units [s], then ωrepresents the angular
frequency and has units [rad/s]. To shorten the notation, we write formally ˆ
f(ω):=
F{f}(ω). On a subspace of L1(R), namely the space S ⊂ L1(R)of so-called Schwartz
functions, the Fourier transformation satisfies a number of important properties. A rigorous
definition of Scan be found e.g. in [Jan71], but it is sufficient for our purpose to think of
it as the space of infinitely differentiable functions decaying rapidly at infinity. For any
function f∈ S it holds that ˆ
f∈ S and the inverse Fourier transform can be defined by
F−1{ˆ
f}(t) = 1
2πZR
ˆ
f(ω)e+iωtdω =f(t).
Moreover, given two functions f, g ∈ S, the following properties of the Fourier transform
follow from its definition:
(2.2) F{f∗g}(ω) = ˆ
f(ω)ˆg(ω)and F{f0}(ω) = +iω ˆ
f(ω),
where f0denotes the first-order derivative of f. The Fourier transform is not restricted
to functions in L1(R)and can be extended to the space of tempered distributions, i.e.,
to the dual space of S. For a rigorous definition of these concepts and for the proof of
properties analogous to (2.2) satisfied by the Fourier transform of tempered distributions,
see e.g. [Jan71].
It is physically meaningful to expect that Eand Hstay bounded at all times since the
system contains no energy source. Therefore, we can assume Eand Hto be tempered
4 M. FROIDEVAUX
distributions and we can apply a Fourier transformation in time to the Maxwell equations
(2.1). This yields to the following constrained eigenvalue problem (EVP):
(2.3) 01
(x,ω)∇×
−1
µ0∇× 0!ˆ
E(x, ω)
ˆ
H(x, ω)= +iω ˆ
E(x, ω)
ˆ
H(x, ω)
∇·((x, ω)ˆ
E(x, ω)) = 0,
∇·(µ0ˆ
H(x, ω)) = 0,
where (x, ω):=0(1 + ˆχE(x, ω)). The first-order (order of derivation in space) EVP
of dimension 6 in (2.3) can be turned into a second-order EVP of dimension 3 simply
by substituing ˆ
E= + 1
iω
1
∇ × ˆ
Hin the second line or, similarly, by substituing ˆ
H=
−1
iω
1
µ0∇× ˆ
Ein the first line. This yields two equivalent formulations of the EVP:
(2.4) ∇×1
(x, ω)∇× ˆ
H(x, ω)=µ0ω2ˆ
H(x, ω)
and
(2.5) ∇×∇× ˆ
E(x, ω) = µ0(x, ω)ω2ˆ
E(x, ω).
Remark that in these two formulations, the constraints ∇ · (ˆ
E)=0and ∇· (µ0ˆ
H)=0
are automatically satisfied since
∇·(∇×f)≡0
for any function ftwice continuously differentiable.
In order to simplify the equations, we restrict our solution space to separable electro-
magnetic fields, i.e. we take the ansatz
(2.6) H(x, t) = H(x)f(t)and E(x, t) = E(x)g(t),
where fand gare two nonzero tempered distributions.
Incorporating the first ansatz into Equation (2.4) and simplifying the terms ˆ
f(ω)ap-
pearing on both sides of the equation yields
(2.7) ∇×1
(x, ω)∇×H(x)=µ0ω2H(x).
Similarly, Equation (2.5) can be rewritten as
(2.8) ∇×∇×E(x) = µ0(x, ω)ω2E(x).
2.2.1. The two-dimensional case. In this subsection, we assume all properties of the system
to be periodic in a two-dimensional (2D) plane and to be constant along the direction
perpendicular to this plane. Without loss of generality, we assume that the plane with
periodic properties is spanned by x1and x2, and that the the direction with constant
properties is along x3. Using the short notation ∂j=∂
∂xjfor j= 1,2,3, the assumption of
constant material properties along x3writes ∂3= 0, and implies ∂3E=∂3H= 0 since no
parameter in Equations (2.7) and (2.8) depends on x3.
Let us introduce some 2D equivalents of the standard three-dimensional (3D) differential
operators, as well as some notation. We consider a 3D vector u= (u1, u2, u3)T, a 2D vector
v= (v1, v2)T, and a scalar function w, and we define
u:=u1
u2, v⊥:=−v2
v1,
grad2Dw:=∂1w
∂2w,div2Dv:=∂1v1+∂2v2,
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 5
curl2Dw:=−(grad2Dw)⊥=∂2w
−∂1w,curl2Dv:=−div2D(v⊥) = ∂1v2−∂2v1,
laplace2Dw:=∂1∂1w+∂2∂2w=−curl2Dcurl2Dw=div2Dgrad2Dw.
It follows that
∇×u=curl2Du3
curl2Du+∂3u⊥
0.
With these tools we can easily show that, in the 2D case, (2.7) becomes
curl2D1
curl2DH1
H2
curl2D1
curl2D(H3)
=µ0ω2H
and (2.8) becomes
curl2Dcurl2DE1
E2
laplace2D(E3)
=µ0ω2E.
Therefore the equations including H3and E3are decoupled from the rest of the equations.
This decoupling can be best expressed in terms of the transverse electric (TE) and trans-
verse magnetic (TM) modes, which we introduce now. Let us start from (2.3) again and
add the effects of the continuous translational symmetry along x3. We obtain
1
curl2D(ˆ
H3)
1
curl2Dˆ
H1
ˆ
H2
1
−µ0curl2D(ˆ
E3)
1
−µ0curl2Dˆ
E1
ˆ
E2
=iω
ˆ
E1
ˆ
E2
(ˆ
E3)
ˆ
H1
ˆ
H2
(ˆ
H3)
where the first and fourth block lines, which define the TE mode (ˆ
ET
1,ˆ
ET
2,0,0,0,ˆ
HT
3)T, are
decoupled from the second and third block lines, which define the TM mode (0,ˆ
HT
1,ˆ
HT
2,
ˆ
ET
3,0,0)T. Combining the equations describing the TE mode yields
(2.9)
(2.10)
−1
µ0
curl2D
1
iω
1
curl2Dˆ
H3=iω ˆ
H3,
1
curl2Dˆ
H3=iω ˆ
E1
ˆ
E2,
and combining the equations describing the TM mode yields
(2.11)
(2.12)
1
curl2D
1
iω
1
−µ0
curl2Dˆ
E3=iω ˆ
E3,
1
−µ0
curl2Dˆ
E3=iω ˆ
H1
ˆ
H2.
It follows that the components ˆ
E1,ˆ
E2,ˆ
H1,ˆ
H2can be easily computed from the solutions
ˆ
H3,ˆ
E3of (2.9) and (2.11).
6 M. FROIDEVAUX
Moreover, the constraints in (2.3) are automatically satisfied in the formulation of Equa-
tions (2.9)–(2.12) since
∇·(ˆ
E) = ∇· 1
iω curl2Dˆ
H3
1
µ0ω2curl2Dcurl2Dˆ
E3!=∇·∇× 1
µ0ω2curl2Dˆ
E3
1
iω ˆ
H3!= 0
and
∇·(µ0ˆ
H) = ∇· 1
−iω curl2Dˆ
E3
curl2D1
ω2curl2Dˆ
H3!=∇·∇×1
ω2curl2Dˆ
H3
−1
iω ˆ
E3= 0.
Finally, note that with the ansatz (2.6), Equations (2.9) and (2.11) can be reformulated
as follows
−div2D(1
grad2DH3) = µ0ω2H3
−laplace2DE3=µ0ω2E3,
or equivalently (for ∂3=∂3H3=∂3E3≡0)
(2.13)
(2.14)
−∇·(1
∇H3) = µ0ω2H3
−∆E3=µ0ω2E3.
The TE mode can be determined by the solution of (2.13) while the TM mode can be
determined by the solution of (2.14).
2.3. Models for the electric permittivity. It has been mentioned in Section 2.1 that we
only consider nonmagnetic materials, i.e. the magnetic permeability in the whole system is
constant and is equal to that of vacuum, namely µ0. The electric permittivity however is
assumed to differ from one material to the other but to stay constant in space within a given
material. Mathematically this means that, if the crystal is composed of Nmat materials
and we denote by Ωjthe domain filled with the j-th material, the electric permittivity of
the crystal can be defined as follows
(x, ω) =
Nmat
X
j=1
j(ω)χΩj(x)
where j(ω)are constant functions in space and χΩjis the indicator function of Ωj.
The functions j(ω)need to be approximated by empirical models. In the remainder
of this subsection, we introduced two of the most commonly used models. Note that we
only consider the dependency of the electric permittivity on the frequency and neglect the
influence of other factors, such as e.g. the temperature.
2.3.1. The frequency-independent model. The simplest model for i(ω)is to assume it to
be constant in frequency, i.e.
i(ω) = (i)
r0
where (i)
rdenotes the relative permittivity of the i-th material and 0the electric permit-
tivity of vacuum.
With this model, the EVPs (2.7), (2.8), (2.13), and (2.14) become linear in the eigenvalue
λ:=ω2.
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 7
2.3.2. The Drude-Lorentz model. We follow here the derivation of the Drude-Lorentz model
presented in [Jac99, p.309-312]. We consider an electron bound to a kernel in a crystal
and submitted to an electric field. We model the interaction between the electron and
the kernel with a spring of stiffness kand of damping constant γ. We make the following
approximations:
•The difference between the local and applied electric fields is neglectable, i.e. the
material has low density;
•Magnetic forces are neglectable;
•The amplitudes of electron oscillations are small enough to evaluate the field at a
constant average position ¯x;
•The applied field Eis assumed to vary harmonically in time.
We denote by q=−ethe charge of the electron and by xits position (with reference
point at the equilibrium point of the spring). Newton’s equation of motion writes in this
case
m¨x=−kx −γ˙x+qE
After defining ω0:=pk/m we obtain
(2.15) m(¨x+γ
m˙x+ω2
0x) = −eE.
Under the assumption that Eis time-harmonic, i.e. with the ansatz E(¯x, t) = E(¯x)eiωt
for a certain ωin C, the Fourier transformation of (2.15) yields to
F{x(t)}(ν) = 2πeE(¯x)
m
1
ω2
0−ν2+iνγ δ(ω−ν) = C(ω)δ(ω−ν)
where C(ω):= 2πeE(¯x)(m(ω2
0−ω2+iωγ))−1is constant in ν. Therefore x(t)is time
harmonic with (angular) frequency ωas well and satifies
x(t) = −e
m
1
(ω2
0−ω2) + iωγ E(¯x, t).
It follows that the dipole moment p:=qx equals
p=e2
m
E
ω2
0−ω2+iωγ .
Let us now consider the dipole moment created by multiple electrons having possibly
different binding frequencies. We denote by Nthe number of molecules per unit volume
and by Zthe number of electrons per molecule. We suppose that there are flelectrons per
molecule with binding frequency ωland damping constant γl, and therefore that Plfl=Z.
We call plthe dipole moment created by the flelectrons corresponding to the binding
frequency ωl. Given the relationships
0
= 1 + ˆχE
and
NX
l
flˆpl=0ˆχEˆ
E
8 M. FROIDEVAUX
we deduce that
0
= 1 + e2N
m0X
l
fl
ω2
l−ω2+iωγl
(2.16)
= 1 + e2N
m0X
l
fl(ω2
l−ω2)
(ω2
l−ω2)2+ω2γ2
l−iωγl
(ω2
l−ω2)2+ω2γ2
l.(2.17)
The damping constant γlare generally small compared to the binding frequencies ωl,
therefore is approximately real for almost all frequencies ω. We notice that the j-th term
in the sum in (2.16) becomes purely imaginary when ω=ωl. Moreover, at this point, the
imaginary part becomes significantly larger in magnitude than at points away from the
resonance.
2.3.3. Physical properties of the electric permittivity. In the previous section, we have seen
that taking into account the damping from the crystal leads to a complex permittivity.
Since the damping is a non-conservative force dissipating kinetic energy into other types
of energy, e.g. heat or sound, a non-real permittivity corresponds to a dissipative system.
Actually, it will be discussed in subection 4.2 that the eigenvalue problems 2.8 and
2.7 are self-adjoint and have their spectrum (i.e. the possible values of ωin the ansatz)
contained in Rif and only if is a real function.
Under the assumption that Evaries harmonically in time, i.e.
E(x, t) = <E(x)eiωt=<E(x)ei(<(ω)+i=(ω))t
=E(x)e−=(ω)t<ei<(ω)t,
the complex part of the frequency describes the decay of the eigenmode over time. In this
setting, in particular with this sign convention in the ansatz, ωcan only be located in the
upper half of the complex plane, since no energy source is present in the crystal.
3. The Floquet transformation & Bloch theorem
3.1. Some fundamentals concepts in crystallography. In this section we introduce
some important concepts in the theory of periodic structures. The reader can refer to the
standard textbook [Kit05] for a more extensive description of these concepts in a physical
context. A detailed presentation of these same concepts in a mathematical context can be
found in [Sim13].
An ideal crystal is formed of two or more components that are arranged in a periodic
manner, i.e. as an infinite repetition of an identical pattern. In this paper we consider
an ideal n-dimensional crystal composed of two different materials, with n=2 or n=3.
Mathematically, one uses the concept of lattice to represent such a periodic structure, and
multiple definitions of this concept can be found in the literature. We choose to use the
following one for its mathematical precision:
Definition 3.1. Let a1, . . . , an∈Rnbe linearly independent. Then the set
G:={g∈Rn|g=
n
X
i=1
νiaifor some νi∈Z}.
is called a (Bravais) lattice and the vectors a1, . . . , anare called the primitive lattice vectors.
With this definition, it follows that any periodic structure can be expressed as a lattice
of repeating objects [Sim13]. This statement is illustrated informally in Figure 3.1.
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 9
Figure 3.1. Decomposition of a periodic structure into a tiling of repeti-
tive units located on a lattice.
Any lattice Gis the result of the tiling of Rnby infinitely many identical tiles called unit
cells, and a unit cell is called primitive if it contains exactly one lattice point. Equivalently,
a primitive unit cell can be thought of as a domain in Rnof smallest volume that can tile
Gusing only translations. In the remainder, we use a particular type of primitive unit cell,
namely the Wigner-Seitz cell defined as follows:
Definition 3.2. Let Gbe a lattice in Rnand g∈G. The Wigner-Seitz cell of Gat the
point gis
Ωg:={r∈Rn| kr−gk<kr−˜gkfor all ˜g∈G, ˜g6=g},
where k·kis the standard Euclidean norm in Rn. One may refer simply to the Wigner-
Seitz Ωof Gwithout mentioning any point gin order to denote tacitly the Wigner-Seitz
cell centered around the origin (0∈G).
To any Bravais lattice Gin Rncorresponds a reciprocal lattice G0defined by
G0:={g0∈Rn|for all g∈G, there exists m∈Zsatisfying g0·g= 2πm}.
It can be shown (see e.g. [Kit05]) that the primitive lattice vectors b1, . . . , bnof G0can
always be chosen in a such way that they satisfy the following constraints: For all j, l =
1, . . . , n,
(3.1) bj·al= 2πδjl,
where δjl denotes the Kronecker delta, and a1, . . . , andenote the primitive lattice vectors
of G.
The Wigner-Seitz cell of G0is called the Brillouin zone of G. In addition to the discrete
translational symmetries represented by the lattice vectors, a lattice Gcan be invariant
under other sorts of symmetries, such as e.g. rotational or axial symmetries. In such a
case, the Brillouin zone reduced by all the symmetries of the lattice is called the irreducible
Brillouin zone Kof G(see Figure 3.2).
3.2. Learning from the electron in a crystal. In this susbsection we recall physical
results historically arising from the analysis of an electron in a perfect crystal.
3.2.1. Translational operators. For any ∆x∈Rnwe denote by T∆x:L2(Rn)→L2(Rn)
the translational operator shifting a function to the right by ∆x, i.e.
T∆xv(x) = v(x−∆x)
for all v(x)∈L2(Rn). Its dual operator with respect to the inner product
(v, w)L2(Rn):=Z
Rn
v(x)w(x)dxfor all v, w ∈L2(Rn)
10 M. FROIDEVAUX
Irreducible
Brillouin
zone
K
of G.
a1
a2
Square lattice G with
lattice vectors a1, a2 and
Wigner-Seitz cell Ω.
Ω
2π/a1
2π/a2
Reciprocal lattice G’
of G with Brillouin
zone Ω’ of G.
Ω’
Ω’
Figure 3.2. Illustrations of some fundamental concepts in crystallography.
is
T∗
∆x:=T−∆x,
since for all v, w, ∈L2(Rn)it holds
(T∗
∆xv, w)=(v, T∆xw) = Z
Rn
v(x)w(x−∆x)dx=Z
Rn
v(˜x+ ∆x)w(˜x)d˜x= (T−∆xv, w).
Thus
T∆xT∗
∆xv(x) = T∆xv(x+ ∆x) = v(x)
for all v∈L2(Rn)and therefore the translational operators is unitary. It follows in
particular that its spectrum is contained on the unit circle in the complex plane.
3.2.2. Bloch theorem.
Lemma 3.3. Given a n-dimensional lattice Gwith lattice vectors {aj}n
j=1, if a function
ψ∈L2(Rd)is an eigenfunction of all translational operators {Taj}n
j=1, then ψis a Bloch
wave, i.e., it is of the form
ψ(x) = eik·xuk(x)
for some periodic function uk(x)with periodicity vectors {aj}n
j=1 and for some vector kin
the Brillouin zone Ω0.
Proof. If ψ(x)is an eigenfunction of all translational operators {Taj}n
j=1, then for all j=
1, . . . , n there exists a real constant θj∈[−1
2,1
2)such that
ψ(x+aj) = e2πiθjψ(x),
since Tajis unitary. Let us denote by {bj}n
j=1 the lattice vectors of the reciprocal lattice
G0satisfying Equation (3.1), and define
k:=
n
X
j
θjbj.
Note that k∈Ω0per construction. Then the function
uk(x):=e−ik·xψ(x).
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 11
satisfies
uk(x+aj) = e−ik·xe−ik·ajψ(x+aj) = e−ik·xe−i2πθjei2πθjψ(x) = e−ik·xψ(x) = uk(x)
for all j= 1, . . . , n, which terminates the proof.
Theorem 3.4 (Bloch theorem).The wavefunction describing the state of electron in a
crystal can be expressed in a basis formed of Bloch waves. Moreover, these Bloch waves are
energy eigenstates.
Proof. Due to the discrete symmetry of the problem, the Hamiltonian and the translation
operators commute with each others. Therefore they can be diagonalized in a common
basis. This basis is thus formed of functions that are energy eigenstates and have the
form of Bloch waves. The expansion of the wavefunction in this basis yields to a linear
combination of Bloch waves, which terminates the proof.
3.3. Floquet transformation. The Bloch theorem is very well-known in solid state
physics. In the mathematics community, this result is more commonly known as a part
of the Floquet theory. Eventhough the Bloch and the Floquet theories are analogous,
their derivations differs a bit. In this subsection we recall the definition of the Floquet
transformation, as well as some of the main results related to it.
Definition 3.5. For any function u∈L2(Rn), we define the Floquet transform of uwith
respect to the lattice Gas the function
F{u}(k, x):=X
g∈G
u(x−g)e−ik·(x−g).
For the sake of brevity, we use the notation uk(x):=F{u}(x, k). Some important
properties follow directly from this definition:
•ukis G−periodic in xand quasi-G0-periodic in k, i.e.
uk(x+g) = uk(x)for all g∈G,
uk+g0(x) = e−ig0·xuk(x)for all g0∈G0,
•given a G-periodic function (x)such that u ∈L2(Rn), it holds
F{u}(k, x) = (x)uk(x),
•for j= 1, . . . , n,
(3.2) F{∂ju}(k, x)=[∂j+ikj]uk(x),
where ∂jdenotes the partial derivative with respect to the j-th component of x
and kjis the j-th component of k.
The following theorem ensures that the Floquet transformation is invertible when the
codomain is properly defined in terms of the irreducible Brillouin zone Kof the lattice G.
Theorem 3.6 ([Kuc01]).Let us denote by Tnthe n-dimensional torus corresponding to the
Wigner-Seitz cell Ωof Gwith periodic boundary conditions. The Floquet transformation
with respect to G, considered as a mapping
F:L2(Rn)→L2(K, L2(Tn))
u→(k→uk(·)),
is isometric and its inverse is given by
u(x) = F−1{uk}(x) = 1
|K| Z
K
uk(x)eik·xdk.
12 M. FROIDEVAUX
3.4. Link between Bloch waves and Floquet transformations. Similarly to the
state of an electron in a crystal, the state of a photon in a crystal, represented by the
electromagnetic field, can be expressed as a linear combination of Bloch waves. Indeed, let
udenote either the electric field Eor the magnetic field H. Then
u(x) = F−1{F{u}}(x) = Z
K
1
|K|eik·xuk(x)dk =Z
K
1
|K|ψk(x)dk,
where ψk(x):=eik·xukxis a Bloch wave.
3.5. The Floquet transformation applied to the Maxwell equations. In the re-
mainder, G⊂Rnrepresents the lattice underlying the spatial distribution of , i.e. (x, ω)
is G-periodic in x. In this subsection, we apply a Floquet transformation with respect to
Gto the Maxwell EVPs for 3D and 2D crystals. In view of the property (3.2) we introduce
the notation ∇k:=∇+ik for any k∈R3and ∇k:= [∂1, ∂2]T+ik for any k∈R2.
3.5.1. 2D photonic crystals. Applying the Floquet transformation to (2.13) and (2.14) for
some kin the irreducible Brillouin zone K ⊂ R2gives
(3.3)
(3.4)
−∇k·(1
∇kH3,k) = µ0ω2H3,k
−∇k·∇kE3,k =µ0ω2E3,k,
respectively, where (3.3) and (3.4) are now two equations on the Wigner-Seitz cell Ω, with
periodic boundary conditions on the functions H3,k(x),E3,k(x)and (x).
3.5.2. 3D photonic crystals. Similarly as in the 2D case, applying the Floquet transforma-
tion to (2.7) and (2.5) for some kin the irreducible Brillouin zone K ⊂ R3gives
∇k×1
(x, ω)∇k×Hk(x)=µ0ω2Hk(x)
∇k×∇k×Ek(x) = µ0(x, ω)ω2Ek(x).
respectively, which are now two equations on the Wigner-Seitz cell Ω, with periodic bound-
ary conditions on the functions Hk(x),Ek(x)and (x).
4. Finite Element Formulation of the EVP
In the remainder we focus our analysis on 2D photonic crystals. Equations (3.3) and (3.4)
are very similar and so are the methods to solve them as well. Therefore, we concentrate
the detailed analysis on (3.4) and keep in my mind that the analysis of (3.3) is analogous.
Moreover, from now on, we consider the easiest empirical model for the electric permittivity,
i.e. we assume that =(x) = Pj=1,2jIΩjis independent of the (angular) frequency ω.
In this case, (3.4) depends only quadratically on ω, so we define the new variable λ:=ω2.
4.1. Weak and discrete formulations. We aim at solving the EVP (3.4) by means of
the finite element method. To do so, we first define the weak formulation of the EVP,
namely: For all k∈ K, find (u, λ)∈V×Csuch that
(4.1) ak(u, v) = λb(u, v)for all v∈V,
where V:=H1
per(Ω) = H1(T2)for T2the 2-dimensional torus corresponding to the Wigner-
Seitz cell Ω, and where a:V×V→Cand b:V×V→Care defined by
ak(u, v):=Z
Ω
∇ku·∇kvdx
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 13
and, respectively,
b(u, v):=Z
Ω
µ0(x)uv dx.
Remark 4.1.We have used the following integration by parts formula for functions in V:
ZΩ−∇·∇kuv dx=ZΩ∇ku·∇vdx−Z∂Ω
(∇ku·ˆn)vds=ZΩ∇ku·∇vdx
where ˆndenotes the the outward unit surface normal to ∂Ωand where R∂Ω(∇ku·ˆn)vds= 0
is to be understood in the sense of traces. This yields
ZΩ−∇k·∇kuv dx=ZΩ∇ku·∇v−ik ·∇kuv dx=ZΩ∇ku·∇kvdx.
With the aim of solving Problem (4.1) numerically, we define a subspace Vh⊂Vof
dimension N < ∞and Nbasis functions {φj}j=1,...,N such that
(4.2) Vh=span{φ1,...φN}.
It follows that any vh∈Vhcan be expressed as a linear combination of the basis functions,
i.e. there exist scalars v1, . . . , vNsuch that vh=
N
P
j=1
vjφj. Thus we can define a discretized
problem corresponding to (4.1) as follows: For all k∈ K, find (λh, uh)∈C×Vhsuch that
(4.3) Aku=λhµ0Bu.
where u:= [u1, . . . , uN]T, and Ak, Bare defined by (Ak)j,i :=ak(φi, φj)and (B)j,i :=
b(φi, φj), for i, j = 1, . . . , N.
Before discussing about methods to solve Problem (4.1) and (4.3), let us analyze the
properties of the linear forms and the corresponding matrices defining the EVPs.
4.2. Properties of the EVPs. Let us start this subsection with some comments and
results concerning the considered vector spaces.
4.2.1. Vector spaces considerations. In the Hilbert spaces H:=L2(Ω) and V=H1(Ω) ⊂
Hwe consider the standard norms
kvkH:=ZΩ|v(x)|2dx
for all v∈Hand
k·kV:=k·kH1(Ω) =k∇·k2
H+k·k2
H1/2.
Remark 4.2.For any u∈V, it holds ∇u∈L2(Ω)2. Therefore the notation k∇ukL2(Ω)
is a bit abusive. In fact, for any w, w0∈L2(Ω)2we use the conventional notation
(w, w0)L2(Ω) :=ZΩ
w(x)·w0(x)dx
where ·denotes the standard dot product, namely w(x)·w0(x):=Pi=1,2wi(x)w0
i(x)for
all w(x), w0(x)∈C2. Additionally,
kwk2
L2(Ω) := (w, w)L2(Ω) =ZΩ|w|2dx
where |w|:=Pi=1,2wiwi1/2denotes the standard norm in C2.
14 M. FROIDEVAUX
When (x)is a real-valued function, we can define another inner product in Hthat will
be more handy for the spectral analysis of the problem.
Lemma 4.3. Let 0, ∞be some positive contants and let (x)be a real-valued function
satisfying 0< 0≤(x)≤∞for all x∈Ω. Then the sesqui-linear form bdefines an
inner product in H. Moreover, the norm k·kbderived from this inner product is equivalent
to k·kHin H.
Proof. The form bpossesses the following properties:
•Positivity: b(u, u) = RΩ(x)|u|2dx≥0kuk2
H≥0and b(u, u)=0if and only if
u≡0,
•Conjugate symmetry: b(u, v) = RΩ(x)uv dx=RΩ(x)vu dx=b(u, v)since
(x)∈R,
•Linearity in the first argument: b(αu, v) = αb(u, v)for all α∈C.
Moreover, it holds that ∞kuk2≥b(u, u)≥0kuk2, which shows the equivalence of the
norms.
In the remainder, we assume (x)to be a real-valued function satisfying 0< 0≤(x)≤
∞for all x∈Ω. Hence, we can define the inclusion map I:V→V∗by means of the
Gelfand triple
V ,→H∼
=H∗,→V∗,
i.e., I=ιH∗→V∗jHιV→Hwhere jH:H→H∗is the Riesz isomorphism with respect to
the inner product b,ιV→His the embedding of Vinto H, and ιH∗→V∗= (ιV→H)∗. An
important property of this particular Gelfand triple is that the imbedding of Vinto His
compact, see e.g. [Eng10] and references therein. It follows that the embedding of H∗into
V∗is also compact [Hac92, Ch. 6].
4.2.2. Definition of the operator equation. In this subsection, we show that Problem (4.1)
can be reformulated as an equivalent problem in the dual space V∗. e start with a lemma
proving Young’s inequality in the specific case of complex vectors.
Lemma 4.4 (Young’s inequality).Consider a,b∈C2. Then, for every δ > 0we have an
estimate of the form
|a·b+a·b| ≤ 1
δ|a|2+δ|b|2.
Proof. For any two vectors c,d∈C2, the following estimates hold:
0≤ |c+d|2= (c+d)·(c+d) = |c|2+|d|2+c·d+c·d,
0≤ |c−d|2= (c−d)·(c−d) = |c|2+|d|2−c·d−c·d.
As a result, we have |c·d+c·d|≤|c|2+|d|2. The claim then follows by setting c=a/√δ,
and d=b√δ.
By means of the Young inequality, we can now show the boundedness of ak:
Lemma 4.5. For all k∈ K, the sesquilinear form ak(·,·)in (4.1) is bounded (or, equiva-
lently, continuous) in V, with constant CM= 2 max(1,|k|2). Moreover, akis Hermitian.
Proof. By means of the Cauchy-Schwarz inequality we obtain for all u, v ∈Vthat
|ak(u, v)| ≤ ZΩ|∇ku(x)||∇kv(x)|dx
≤ZΩ|∇ku(x)|2dx1/2ZΩ|∇kv(x)|2dx1/2
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 15
Lemma 4.4 with δ= 1 yields
ZΩ|∇ku(x)|2dx≤ZΩ|∇u(x)|2+|ku(x)|2+iku∇u(x) + iku∇u(x)
2dx
≤ k∇uk2
H+|k|2kuk2
H+kkuk2
H+k∇uk2
H.
It follows that
|ak(u, v)| ≤ CMkukVkvkV
where CM= 2 max(1,|k|2). Moreover, it is straightforward to see that ak(u, v) = ak(v, u),
which implies the Hermiticity of ak.
Hence, we can define an operator Akmapping from the space Vto its dual V∗by the
equality
hAku, ·iV∗,V :=ak(u, ·)
for all u∈V. The EVP (4.1) can then be reformulated as an equation in V∗:
(4.4) Akuk=λIuk.
4.2.3. Properties of the resolvent. Instead of the eigenvalue problem (4.4) it is more con-
venient for the analysis to look at the equivalent shifted problem:
(Ak+σI)uk= (λ+σ)Iuk
where σ > 0is a fixed parameter. In this section we investigate the properties of the
resolvant (Ak+σI)−1of Ak. To this end, we demonstrate that the form ak,σ :=ak+σb
corresponding to the operator Ak,σ :=Ak+σIis an elliptic sesquilinear form.
Lemma 4.6. For any σ > 0, the sesquilinear form ak,σ is elliptic with constant
Cm(σ) = (min{1, σ0},if k= 0,
min nσε0
2|k|2+σε0,σ
2ε0o,if k6= 0.
Proof. If k= 0, then a0,σ(u, u) = RΩ|∇u|2+σ(x)|u|2dx≥min{1, σ0}kuk2
V. If k6= 0, we
observe that ak,σ(u, u) = RΩ|(∇+ik)u|2+σ(x)|u|2dx∈R+. Moreover,
ak,σ(u, u)≥k∇uk2
H+ (σ0+|k|2)kuk2
H−ZΩiku ·∇u+iku ·∇udx
and thanks to Lemma 4.4 we obtain :
ZΩiku ·∇u+iku ·∇udx≤|k|2kuk2
H
δ+δk∇uk2
H
for any δ > 0. Therefore,
ak,σ(u, u)≥(1 −δ)k∇uk2
H+σε0+|k|2−|k|2
δkuk2
H,
and choosing δ=|k|2
|k|2+σ
2ε0
<1yields
ak,σ(u, u) = (1 −|k|2
|k|2+σ
2ε0
)
| {z }
>0
k∇uk2
H+σ
2ε0
|{z}
>0
kuk2
H≥Cm(σ)kuk2
V
with Cm(σ) = min nσε0
2|k|2+σε0,σ
2ε0o.
16 M. FROIDEVAUX
Corollary 4.7. For any σ > 0, the operator Ak,σ = (Ak+σI)is invertible and the
inverse is bounded with constant Cm(σ)−1, i.e. Ak,σ−1∈ L(V∗, V )satisfies
kAk,σ−1kV∗→V:= sup
w∈V∗
kwkV∗=1
kAk,σ−1wkV≤1
Cm(σ).
Proof. This follows directly from the Lax-Milgram lemma, see e.g. [Hac92, Lem. 6.98].
Corollary 4.8. For any σ > 0, the sesquilinear form ak,σ defines an inner product in V
and the corresponding norm kvkak+σb:=pak(v, v) + σb(v, v)is equivalent to the standard
norm kvkV.
Proof. The positivity of ak,σ follows from Lemma 4.6, and the sesquilinearity as well as
the conjugate symmetry are easily verified. Moreover Lemma 4.5 yields
kvkak+σb≤qCMkvk2
V+σb(v, v)≤qCMkvk2
V+σ∞kvk2
H≤pCM+σ∞kvkV
and by Lemma 4.6
kvkak+σb≥pCm(σ)kvkV
for all v∈V.
Many important properties of the spectrum of Akcan now be proven by applying the
Riesz-Schauder therory on the resolvant Rσ:V→V,R=A−1
k,σI. These are gathered in
the following lemma.
Lemma 4.9. The spectrum of Akconsists in at most countably many real eigenvalues
with finite geometrical multiplicities, which can only accumulate at infinity. Moreover, all
eigenvalues are nonnegative.
Proof. From the Hermiticity of akfollows directly that Akis Hemitian. Therefore all
eigenvalues λof Akmust be real since
λ=hAku, uiV∗,V
hIu, uiV∗,V
=hAku, uiV∗,V
hIu, uiV∗,V
=λ,
where udenotes the eigenvector corresponding to λ. The rest of the proof follows from
[Hac92, Satz 6.108]. In particular,it is shown there that for all λ∈Cone of the following
alternatives holds:
(1) (A−λI)−1∈ L(V∗, V ),
(2) λis an eigenvalue.
In Corollary 4.7 it was shown that (A−λI)−1exists and is bounded for all λ=−σ, where
σ > 0. Hence, all eigenvalues must be nonnegative.
Corollary 4.10. Let Akand Bbe the matrices in the discretized problem (4.3). Then,
Akis Hermitian positive semi-definite and B, as well as (Ak+B), are Hermitian positive
definite.
Proof. It is straightforward to see that (Ak)j,i =ak(φi, φj) = ak(φj, φi) = (Ak)i,j and
similarly for B. Moreover, by Lemma 4.6 it holds for all σ > 0and for all v= [v1, . . . , vN]
with vh=PN
ivi
v∗Akv=a(vh, vh)≥Cm(σ)kvhkV−σb(vh, vh)→0as σ→0,
where the superscript ∗denotes the conjugate transpose of a vector or a matrix Since this
inequality holds for σ > 0arbitrarily small, to follows that v∗Akv≥0, and thus Akis
positive semi-definite. The rest of the proof follows easily.
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 17
5. Error estimators
We wish to approximate the solution of the discretize problem (4.3) by using an iterative
method. In the remainder, udenotes the exact eigenfunction of (4.1), uhdenotes the
exact eigenfunction of the discretized problem (4.3), and euh=PN
i=1 euiφidenotes the
approximated eigenfunction computed by msteps of an iterative method (e.g. Lanczos
method or LOBPCG [GVL96, Kny01]).
In order to be able to evaluate the quality of the approximate solution euhwe need
to construct an estimator for the error ku−euhkV. We start by applying the triangular
inequality
ku−euhkV=ku−uh+uh−euhkV≤ ku−uhkV+kuh−euhkV
and observe that the error can be separated into the discretization error ku−uhkVand
the algebraic error kuh−euhkV.
An estimator for the error introduced by the discretization of the linear EVP (3.4) to
obtain (4.3) has been developed in [GG12]. In that paper it is shown that an error estimator
ηdis, which can be computed a posteriori, gives a reliable and efficient approxmation of the
error ku−euhkak+σb where σ= 1. Since the norms k·kVand k·kak+σb are equivalent
(see Lemma 4.8), ηdis is also a reliable and efficient estimator of the error measured in the
standard norm in V. However, the equivalence constants Cm(σ)and CMmay be much
bigger or much smaller than 1, depending on the values of σand k.
In order to enable the error balancing between the discretization and the algebraic errors,
we need an error estimator which can be compared to ηdis in a meaningful manner, i.e., an
estimator which measures the algebraic error in the norm defined by the sesquilinear form
ak,σ on V. This is the subject of the next subsection.
5.1. Algebraic error. We aim at finding a bound for the algebraic errors kuh−˜uhkV
and kuh−˜uhkak+σb, which depend on the convergence of an iterative numerical method.
The convergence analysis of such methods is typically based on the concepts of angle
between vectors. While the definition of the angle between real vectors with respect to
a scalar product is clear, there exist multiple nonequivalent concepts of angle between
complex vectors. Each one of these angles has different interesting properties. We start
by introducing the Euclidean angle between complex vectors, which satisfies the law of
cosines:
Definition 5.1. Let vand wbe two vectors in a complex Hilbert space Vand let (·,·)I
as well as k·kIdenote an inner product in Vand its corresponding norm, respectively.
The (real-valued) Euclidean angle between vand wwith respect to the I-inner-product is
defined as the quantity ∠(E)
I(v, w)∈[0, π]satisfying
∠(E)
I(v, w) = acos <(v, w)I
kvkIkwkI
.
The Euclidean angle indeed satisfies the law of cosines since for all v, w ∈Vit holds
kv−wk2
I= (v−w, v −w)I=kvk2
I+kwk2
I−(v, w)I−(w, v)I
=kvk2
I+kwk2
I−2<(v, w)I=kvk2
I+kwk2
I−2 cos ∠(E)
I(v, w)kvkIkwkI.
The Euclidean angle inherits other properties from the real case, e.g. the angle between
vectors pointing in opposite directions is πsince ∠(E)
I(v, −v) = acos(−1) = π. However,
the Euclidean angle is not invariant under scaling since, for α∈C,∠(E)
I(αv, w)is generally
not equal to ∠(E)
I(v, w).
18 M. FROIDEVAUX
On the other hand, one can define another angle between complex vectors which is
invariant under scaling:
Definition 5.2. Let v, w ∈V,(·,·)Iand k·kIbe as in Definition 5.1. The Hermitian
angle between vand wwith respect to the I-inner-product is defined as the quantity
∠(H)
I(v, w)∈[0,π
2]satisfying:
∠(H)
I(v, w) = acos |(v, w)I|
kvkIkwkI
.
The link between both angles can be expressed as follows:
(5.1)
inf
vθ∈span{v}
wθ∈span{w}
∠(E)
I(vθ, wθ) = inf
θ1∈[0,2π)
θ2∈[0,2π)
∠(E)
I(eiθ1v, eiθ2w) = acos sup
θ∈[0,2π)
<eiθ(v, w)I
kvkIkwkI
= acos 1
kvkIkwkI<(w, v)I
|(v, w)I|(v, w)I
= acos |(v, w)I|
kvkIkwkI
=∠(H)
I(v, w).
Therefore, the Hermitian angle between vand wcan be thought of as the Euclidean angle
between the subspaces spanned by vand w. Moreover, the Hermitian and Euclidean angles
between two vectors v, w are equal if the vectors have the same orientation in V, i.e. if
they satisfy the condition
(v, w)I∈R+.
In the remainder, we consider the algebraic error kuh−˜uhkIwhere the inner prod-
uct (·,·)Irepresents either the standard inner product (·,·)Vin Vor the inner product
(·,·)ak+σblinked to the weak form (4.1). Moreover, both vectors uhand ˜uhare supposed
to be normalized and have the same orientation in V(with respect to the considered inner
product), i.e.
kuhkI=k˜uhkI= 1 and (uh,˜uh)I∈R+.
Under these conditions, the law of cosines combined with the following trigonometric iden-
tity
cos2(φ
2) = eiφ
2+e−iφ
2
2!2
=eiφ +e−iφ + 2
4=1
2(cos(φ) + 1),for all φ∈[0,2π),
yields
kuh−˜uhk2
I= 2 −2 cos ∠(E)
I(uh,˜uh)=41−1
2cos ∠(E)
I(v, w)+1
= 41−cos2∠(E)
I(v, w)
2= 4 sin2∠(E)
I(v, w)
2= 4 sin2∠(H)
I(v, w)
2.
Therefore, we can bound the algebraic error by a term proportionnal to the sinus of the
Hermitian angle:
(5.2) kuh−˜uhkI= 2 sin ∠(H)
I(uh,˜uh)
2≤√2 sin ∠(H)
I(uh,˜uh).
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 19
For small angles, the following expression based on Taylor expansions gives a more accurate
estimation of the error:
kuh−˜uhkI= 2 sin ∠(H)
I(uh,˜uh)
2= 2
∞
X
j=0
(−1)j
(2j+ 1)!∠(H)
I(uh,˜uh)
22j+1
=
∞
X
j=0
(−1)j
(2j+ 1)!∠(H)
I(uh,˜uh)2j+1 +
∞
X
j=1
(−1)j
(2j+ 1)!(1 −1
22j)∠(H)
I(uh,˜uh)2j+1
= sin ∠(H)
I(uh,˜uh) + O(∠(H)
I(uh,˜uh)3)≈sin ∠(H)
I(uh,˜uh).
5.1.1. Equivalence of continuous and discrete norms. In this subsection we establish equiv-
alences between the norms of functions uh, vhin the finite dimensional space Vhdefined
in (4.2) and the corresponding vectors u,v∈Cnexpressed in the basis formed by the
functions {φi}N
i=1. In the remainder of this paper, we assume σ > 0.
For any uh, vh∈Vhit holds
(uh, vh)V=
N
X
i=1
N
X
j=1
vjZΩ∇φi·∇φj+φiφjdx ui=v∗ΦVu
and
kuhk2
V=uhΦVu=:kuk2
ΦV,
where (ΦV)i,j := (φj, φi)V. Moreover,
(uh, vh)b=
N
X
i=1
N
X
j=1
vjZΩ
(x)φiφjdx ui=v∗Bu,
kuhk2
b=u∗Bu=kuk2
B
and
(uh, vh)2
ak+σb=
N
X
i=1
N
X
j=1
vjZΩ∇kφi·∇kφj+σφiφjdx ui=v∗(Ak+σB)u,
kuhk2
ak+σb=u∗(Ak+σB)u=kuk2
(Ak+σB).
Similarly, a discrete norm that is equivalent to a norm in V∗can be found.
Lemma 5.3. For any r∈V∗such that r=Irhfor some rh=
N
P
j=1
rjφj∈Vh, it holds
krkV∗=kBrk(Ak+σB)−1=krkB(Ak+σB)−1B,
where
kfkV∗:= sup
v∈V,
kvkak+σb=1
|f(v)|
for any f∈V∗, and where r= [r1, . . . , rN].
20 M. FROIDEVAUX
Proof. For any v∈Vwe denote by πh(v)the orthogonal projection of vonto Vhwith
orthogonality in the inner product (·,·)b. Then it holds
krkV∗= sup
v∈V,
kvk(ak+σb)=1
|(rh, v)b|= sup
v∈V,
kvk(ak+σb)=1
|(rh, v −πh(v))b+ (rh, πh(v))b|
= sup
v∈V,
kvk(ak+σb)=1
|(rh, πh(v))b|= sup
vh∈Vh,
kvhk(ak+σb)=1
|(rh, vh)b|
= sup
v∈CN,
kvk(Ak+σB)=1
|r∗Bv|= sup
v∈CN,
kvk(Ak+σB)=1
|r∗B(Ak+σB)−1(Ak+σB)v|.
For the sake of brevity, we introduce the notation Ak,σ :=Ak+σBand b:=A−1
k,σBr.
Using the Cauchy-Schwarz inequality [Bre11, Ch. 5.1] yields
krkV∗= sup
v∈CN,
kvkAk,σ =1
|b∗Ak,σv| ≤ sup
v∈CN,
kvkAk,σ =1
kbkAk,σ kvkAk,σ = sup
v∈CN,
kvkAk,σ =1
kbkAk,σ =kbkAk,σ .
The upper bound is attained by setting v=b
kbkAk,σ
which satisfies the constraint kvkAk,σ =
1. Therefore
krkV∗=kbkAk,σ =kA−1
k,σBrkAk,σ =kBrkA−1
k,σ =krkB(Ak+σB)−1B.
We show that the trigonometric functions of angles between vectors in Vhcan also be
expressed in terms of the corresponding discrete vectors in CN.
Theorem 5.4. Let vh, wh∈Vh⊂Vwith vh=
N
P
i=1
viφi,wh=
N
P
i=1
wiφiand v=
[v1, . . . , vN]T,w= [w1, . . . , wN]T. Then
cos ∠(H)
V(vh, wh) = cos ∠(H)
ΦV(v,w)
and
cos ∠(H)
ak+σb(vh, wh) = cos ∠(H)
Ak+σB(v,w)
where
∠(H)
M(v,w):= acos |w∗Mv|
√v∗Mv√w∗Mw∈h0,π
2i
for any matrix M.
Proof.
cos ∠(H)
V(vh, wh) = |(vh, wh)V|
kvhkVkwhkV
=|P
i,j
wi(φj, φi)Vvj|
rP
i,j
vi(φj, φi)VvjrP
i,j
wi(φj, φi)Vwj
=|w∗ΦVv|
√v∗ΦVv√w∗ΦVw
and similarly for the inner product ak+σbon V.
Corollary 5.5. Let vh, wh,v,was in Theorem 5.4. Then
sin ∠(H)
V(vh, wh) = sin ∠(H)
ΦV(v,w)
and
sin ∠(H)
Ak+σB(vh, wh) = sin ∠(H)
ak+σb(v,w).
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 21
5.1.2. Sine theorems. In order to obtain a useful bound on the algebraic error, we need a
computable bound on sin ∠(H)
ΦV(v,w)or sin ∠(H)
Ak+σB(v,w). First, we recall a result for the
(non-generalized) linear eigenvalue problem.
Theorem 5.6 (The sine theorem).Let M∈CN×Nbe a Hermitian positive definite matrix
defining an inner product (v1,v2)Msuch that
(v1,v2)M:=v1∗Mv2
for all v1,v1∈CN. Let A∈CN×Nbe a Hermitian matrix with respect to this inner
product. We denote by {λi}N
i=1 the eigenvalues of Aand by {u(i)}N
i=1 the corresponding
normalized eigenvectors of Ain the norm k·kMinduced by (·,·)M. Given any vector
v∈CNand any scalar e
λ∈C, we order the eigenpairs (λi,u(i))such that |λ1−e
λ| ≤
|λ2−e
λ| ≤ . . . ≤ |λN−e
λ|and we define the residual r:=Ae
v−e
λe
vwhere e
v:=αvis the
scaled vector which is normalized and lies in the same directions as u(1), i.e.
ke
vkM= 1,(e
v,u(1))M∈R+.
Then it holds
sin ∠(H)
M(e
v,u(1))≤krkM
δ.
where δ:=|λ2−e
λ|.
Proof. We follow the proof of Theorem 11.7.1 in [Par98] which proves the result for real
symmetric matrices and we extend it to the case of complex Hermitian matrices.
First recall that the eigenvectors of Aform an orthogonal basis of Cnsince Ais Hermitian
(see the spectral theorem in [Par98, Ch. 1.4]). Therefore we can expand
e
v=
N
X
j=1
(e
v,u(j))Mu(j),
which implies
r= (A−e
λ)(e
v,u(1))Mu(1) + (A−e
λ)
N
X
j=2
(e
v,u(j))Mu(j)
= (λ1−e
λ)(e
v,u(1))Mu(1) +
N
X
j=2
(λj−e
λ)(e
v,u(j))Mu(j).
Since the eigenvectors are orthogonal to each other in the inner product (·,·)M, we can
apply the Pythagoreus theorem, which yields:
krk2
M=|λ1−e
λ|2|(e
v,u(1))M|2+
N
X
j=2 |λj−e
λ|2|(e
v,u(j))M|2
≥ |λ1−e
λ|2cos2∠(H)
M(e
v,u(1)) + |λ2−e
λ|2
N
X
j=2 |(e
v,u(j))M|2
=|λ1−e
λ|2cos2∠(H)
M(e
v,u(1)) + |λ2−e
λ|2(ke
vk2−|(e
v,u(1))M|2)
=|λ1−e
λ|2cos2∠(H)
M(e
v,u(1)) + |λ2−e
λ|2sin2∠(H)
M(e
v,u(1))
≥δ2sin2∠(H)
M(e
v,u(1)).
22 M. FROIDEVAUX
In order to apply Theorem 5.6 to Problem (4.3) we need to reformulate the latter as a
standard (non-generalized) EVP. Since the matrix Bin (4.3) is Hermitian positive definite,
it is also invertible and the EVP can be equivalently written as
B−1
Aku=λu.
Corollary 5.7. Let Ak, Bbe the Hermitian matrices defining the EVP (4.3). We de-
note by {λi}N
i=1 the eigenvalues of the pencil (Ak, B)and by {u(i)}N
i=1 the correspond-
ing normalized eigenvectors of (Ak, B)in the norm k·kAk+σBwhere σ > 0. Given
any vector v∈CNand any scalar e
λ∈C, we order the eigenpairs (λi,u(i))such that
|λ1−e
λ|≤|λ2−e
λ| ≤ . . . ≤ |λN−e
λ|and we define the residual r:=Ake
v−e
λBe
vwhere
e
v:=αvis the scaled vector which is normalized and lies in the same directions as u(1),
i.e.
ke
vkAk+σB= 1,(e
v,u(1))Ak+σB∈R+.
Then it holds
(5.3) sin ∠(H)
Ak+σB(e
v,u(1))≤kB−1
rkAk+σB
δ=krkB−1
(Ak+σB)B−1
δ.
where δ:=|λ2−e
λ|.
Proof. Observe that the matrix B−1
Akis Hermitian in the inner product defined by Ak+
σBsince for any v1,v2∈CN
(B−1
Akv1,v2)Ak+σB=v∗
1AkB−1
(Ak+σB)v2=v∗
1AkB−1
Akv2+σv∗
1Akv2
=v∗
1(Ak+σB)B−1
Akv2= (v1, B−1
Akv2)(Ak+σB).
Therefore Theorem 5.6 can be applied with A:=B−1
Ak,M:= (Ak+σB).
Alternatively, one can follow [HL06, Section 2] and solve the EVP (4.3) with a shift-
invert method. Recall that for any real positive shift σ > 0, the matrix (Ak+σB)is
invertible, see Corollary 4.7. Therefore (4.3) is equivalent to
(5.4) (Ak+σB)−1Bu=νu,
where the eigenvalues of both problems are related via the equation
ν=1
λ+σ.
The matrix (Ak+σB)−1Bis Hermitian in the (Ak+σB)-inner product since for all
v1,v2∈CNit holds
(Ak+σB)−1Bv1,v2Ak+σB=v∗
1B(Ak+σB)−1(Ak+σB)v2=v∗
1Bv2
=v1,(Ak+σB)−1Bv2Ak+σB
We can now apply Theorem 5.6 to (5.4).
Corollary 5.8. Let Ak, Bbe the Hermitian matrices defining the EVP (4.3) and σ > 0
be a real and positive scalar shift. We denote by {λi}N
i=1 the eigenvalues of the Hermitian
pencil (Ak, B)and by {u(i)}N
i=1 the corresponding normalized eigenvectors of (Ak, B).
Given any vector v∈CNand any scalar e
λ∈C, we order the eigenpairs (λi,u(i))such that
1
λ1+σ−1
e
λ+σ≤
1
λ2+σ−1
e
λ+σ≤. . . ≤
1
λN+σ−1
e
λ+σ,
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 23
and we define the residual r:=Ake
v−e
λBe
v, where e
v:=αvis normalized and lies in the
same directions as u(1), i.e.
ke
vkAk+σB= 1,(e
v,u(1))Ak+σB∈R+.
Then it holds
(5.5) sin ∠(H)
Ak+σB(u(1),e
v)≤|λ2+σ|
|λ2−e
λ|krk(Ak+σB)−1.
Proof. Applying Theorem 5.6 with A:= (Ak+σB)−1B,M=Ak+σBand eν:= (e
λ+σ)−1
yields
sin ∠(H)
Ak+σB(u(1),e
v)≤k(Ak+σB)−1Be
v−eνe
vk(Ak+σB)
|eν−ν2|
=|eν|
|eν−ν2|keν−1Be
v−(Ak+σB)e
vk(Ak+σB)−1
=|λ2+σ|
|λ2−e
λ|ke
λBe
v−Ake
vk(Ak+σB)−1.
Remark 5.9.Observe that the norm of the residual rin the upper bound (5.5) can be
rewritten as
krk(Ak+σB)−1=ke
λBe
u−Ake
uk(Ak+σB)−1=ke
λe
u−B−1
Ake
ukB(Ak+σB)−1B.
Therefore, according to Lemma 5.3, krk(Ak+σB)−1can be seen as the norm of
(5.6) res(e
λ, e
u):=I N
X
i=1 B−1
riφi!
in the dual space V∗, i.e.
krk(Ak+σB)−1=kres(e
λ, e
u)kV∗
Note that (5.6) indeed defines a residual, since for all vh=PN
i=1 viφi∈Vhit holds that
hres(e
λ, e
u), vhiV∗,V =ZΩ
(x)
N
X
i=1 B−1
riφivhdx=
N
X
i,j=1 B−1
rivjZΩ
(x)φiφjdx
=
N
X
i,j=1 B−1
rivj(B)j,i =v∗BB−1
r=v∗r=v∗(Ake
u−e
λBe
u)
=ak(euh, vh)−e
λb(euh, vh).
where euh:=PN
j=1 e
ujφj
Remark 5.10.Since (Ak+σB)is Hermitian positive definite we can factorize it with help
of the Choleski decomposition (Ak+σB) = LL∗where Lis invertible. Therefore the term
appearing on the right side of (5.3) can be rewritten as
kB−1
Ake
v−e
λe
vkAk+σB=kL∗B−1
AkL−∗ e
w−e
λe
wk
with L∗e
v=e
w. Observe that the matrix L∗B−1
AkL−∗ is Hermitian since
(L∗B−1
AkL−∗)∗=L−1AkB−1
L=L−1AkB−1
(Ak+σB)L−∗L−1L
=L−1(Ak+σB)B−1
AkL−∗ =L∗B−1
AkL−∗.
24 M. FROIDEVAUX
Therefore, if e
wis the approximation of the eigenvector corresponding to the smallest
eigenvalue of L∗B−1
AkL−∗ returned after msteps of an inverse Lanczos or LOBPCG
method, then kL∗B−1
AkL−∗ e
w−e
λe
wkis directly available as a by-product of each iteration,
see e.g. [GVL96, Thm. 9.1.2] and [Kny01].
Similarly the term on the right hand side of (5.5) can be rewritten as
ke
λBe
v−Ake
vk(Ak+σB)−1=kL−1AkL−∗ e
w−e
λL−1BL−∗ e
wk
and can be obtained as a by-product of a LOBPCG iteration on the Hermitian pencil
(L−1AkL−∗, L−1BL−∗).
Remark 5.11.In Lanczos’ method and LOBPCG, the eigenvalue converges faster than the
eigenvector, therefore it is reasonable to approximate λ≈e
λand λ2≈e
λ2to compute the
error bound.
All results of this section can be summarized in the following theorem:
Theorem 5.12. Let (Ak, B)be the matrix pencil defining the discretized problem (4.3)
and let ube the eigenvector corresponding to the smallest eigenvalue of (Ak, B). We write
Lthe triangular matrix satisfying the Choleski decompostion LL∗=Ak+σBfor some
shift σ > 0. We denote by e
w=:L∗e
uthe approximation of L∗uobtained after miterations
of the LOBPCG method applied to the pencil (L−1AkL−∗, L−1BL−∗), and assume uand
e
uto be scaled in such a way that
kukAk+σB=ke
ukAk+σB= 1,(e
u,u)Ak+σB∈R+.
Let e
λand e
λ2be the approximations of the eigenvalue λclosest to σand, respectively, of the
eigenvalue λ2second closest t σ, of Problem (4.3) obtained from the LOBPCG iterations.
Then the distance between the functions uh=PN
j=1 ujφj∈Vhand euh=PN
j=1 e
ujφj∈Vh
can be approximately bounded by
kuh−˜uhkV≤pCm(σ)−1kuh−˜uhkak+σb=pCm(σ)−1sin ∠(H)
ak+σb(uh,euh) + O(∠(H)
ak+σb(uh,euh)3)
≈pCm(σ)−1sin ∠(H)
ak+σb(uh,euh)≤pCm(σ)−1|λ2+σ|
|λ2−e
λ|kAke
u−e
λBe
uk(Ak+σB)−1
≈pCm(σ)−1|e
λ2+σ|
|e
λ2−e
λ|kL−1AkL−∗ e
w−e
λL−1BL−∗ e
wk
when ∠(H)
ak+σb(uh,euh)is a small angle.
5.2. Numerical experiments. In this section, we compare the error estimator
(5.7) ηalg :=|e
λ2+σ|
|e
λ2−e
λ|kL−1AkL−∗ e
w−e
λL−1BL−∗ e
wk
to the error kweig −e
wk ≈ kuh−˜uhkak+σbwhere weig is the approximated eigenvector
returned by the MATLAB function eig and is therefore assumed here to be a very good
approximation of the exact eigenvector w=L∗u. We observe that ηalg is composed of a
factor ρ=|e
λ2+σ|/|e
λ2−e
λ|depending only on the approximations of the eigenvalues, and
of a factor corresponding to the norm or the residual. In this section, we denote with a
superscript (m)a quantity computed with the approximations e
λ, e
λ2,e
wobtained after m
iterations.
We have tested the error estimator (5.7) on a problem of the form (4.3) with (x, ω)as in
subsection 2.3.1, Ωbeing a square cell of length lΩ= 500·10−9m, Ω1⊂Ωbeing a centered
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 25
disc of radius 0.2lΩ,1(ω)≡0= 8.854187817·10−6µF/m, 2≡8.90. We have set σ= 1,
and Figures 5.1 and 5.2 correspond to k= [0.6888,6.2832]T·106m−1. The problem was
discretized via the C++ library Concepts [FL02] with quadratic Lagrange basis functions.
The dimension of the resulting matrices is 896×896. We solved the Hermitian generalized
EVP with the LOBPCG algorithm.
On Figure 5.1 we see that the error estimator does not give reliable estimation of the
error during the first dozen of iterations, i.e. when the iterative algorithm is far from con-
vergence. However, after approximately 25 iterations, the error estimator starts mimicking
the decrease of the error. We call this behavior asymptotic, as opposed to the transient
behavior happening in the first steps of the iterative process. We observe that there exists
a constant shift Mbetween the graph of the error estimator and the one of the error in
the asymptotic regime, which corresponds to approximately M=13 iterations.
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
0
0.2
0.4
0.6
0.8
1
1.2
1.4
m [nb of iterations]
||weig −e
w||
sin ∠(weig,e
w)
√2 sin ∠(weig,e
w)
ηalg =ρ||res||
Figure 5.1. Comparison of the error with the error estimator
On Figure 5.2 we see that the error, the error estimator, and the norm of the residual
indeed decrease at the same speed, and that the beginning of the asymptotic regime cor-
relates with the stabilization of the factor ρat m≈. This correlation is made clearer in
Figure 5.3, where the rates of change of the factors ρand kreskare plotted. Indeed, we see
that the biggest mat which r(m)
ρ:=|ρ(m)−ρ(m−1)|
ρ(m)is bigger than r(m)
res :=|kres(m)k−kres(m−1)k|
kres(m)k
is also m≈20. In the range m > 20,r(m)
ρbecomes significantly smaller than r(m)
res .
Note that r(m)
ρand r(m)
res can be cheaply evaluated at each iteration step as a complement
of the error estimator ηalg. Comparing both rates of change gives a way to decide whether
the error estimator is reliable or not. The efficiency of the error estimator, however, depends
on the shift M, and it is for the moment not clear how to estimate this quantity.
We have made a qualitatively similar observations for other values of kon the boundary
of the irreducible Brillouin zone. However, the shift Mobserved between the error and the
error estimator in the asymptotic domain does not stay constant with kbut varies within
the range 6−20 iterations.
26 M. FROIDEVAUX
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150
10−20
10−15
10−10
10−5
100
105
m [nb of iterations]
||weig −e
w||
||res||
ρ=|e
λ2+σ|/|e
λ2+e
λ1|
ηalg =ρ||res||
Figure 5.2. Evolution of the factors composing the error estimator
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150
10−20
10−15
10−10
10−5
100
105
m [nb of iterations]
|ρ(m)−ρ(m−1)|/|ρ(m)|
abs(||res(m)||−||res(m−1)||)/||res(m)||
Figure 5.3. Rates of change of the factors composing the error estimator
6. Conclusion
We have developed a residual-based estimator which approximates the algebraic error
resulting from the solution of the discrete generalized linear EVP by an iterative method.
The error estimator approximates the error in the norm arising from the sesquilinear forms
defining the EVP, and can therefore by combined in an error balancing procedure with
already existing estimators of the discretization error. We have shown how to heuristically
determine whether enough iterations have been computed for the error estimator to be
reliable.
ON AFEM FOR THE SIMULATION OF 2D PHOTONIC CRYSTALS 27
References
[BCG06] D. Boffi, M. Conforti, and L. Gastaldi. Modified edge finite elements for photonic crystals.
Numer. Math., 105(2):249–266, 2006.
[Bre11] H. Brezis. Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer,
New York, 2011.
[CDM+17] E. Cancès, G. Dusson, Y. Maday, B. Stamm, and M. Vohralí k. Guaranteed and robust a
posteriori bounds for Laplace eigenvalues and eigenvectors: conforming approximations. SIAM
J. Numer. Anal., 55(5):2228–2254, 2017.
[EKE12] C. Effenberger, D. Kressner, and C. Engström. Linearization techniques for band structure
calculations in absorbing photonic crystals. Int. J. Numer. Meth. Engn, 89(2):180–191, 2012.
[ELT17] Christian Engström, Heinz Langer, and Christiane Tretter. Rational eigenvalue problems and
applications to photonic crystals. J. Math. Anal. Appl., 445(1):240–279, 2017.
[Eng10] C. Engström. On the spectrum of a holomorphic operator-valued function with applications to
absorptive photonic crystals. Math. Models Methods Appl. Sci., 20(8):1319–1341, 2010.
[FL02] P. Frauenfelder and C. Lage. Concepts - An Object-Oriented Software Package for Partial
Differential Equations. ESAIM: M2AN, 36(5):937–951, 2002.
[GG12] S. Giani and I. G. Graham. Adaptive finite element methods for computing band gaps in
photonic crystals. Numer. Math., 121(1):31–64, 2012.
[GVL96] G. H. Golub and C. F. Van Loan. Matrix computations. Johns Hopkins Studies in the Mathe-
matical Sciences. Johns Hopkins University Press, Baltimore, MD, third edition, 1996.
[Hac92] W. Hackbusch. Elliptic Differential Equations: Theory and Numerical Treatment. Springer,
Berlin, 1992.
[HL06] U. L. Hetmaniuk and R. B. Lehoucq. Uniform accuracy of eigenpairs from a shift-invert Lanczos
method. SIAM J. Matrix Anal. Appl., 28(4):927–948, 2006.
[HLM16] T.-M. Huang, W.-W. Lin, and V. Mehrmann. A Newton-type method with non-equivalence
deflation for nonlinear eigenvalue problems arising in photonic crystal modeling. SIAM J. Sci.
Comput., 38(2):B191–B218, 2016.
[Jac99] J. D. Jackson. Classical electrodynamics. John Wiley & Sons, Inc., New York-London-Sydney,
third edition, 1999.
[Jan71] L. Jantscher. Distributionen. Walter de Gruyter & Co, Berlin, 1971.
[Kit05] C. Kittel. Introduction to Solid State Physics. John Wiley & Sons, Inc., Hoboken, NJ, eighth
edition, 2005.
[Kny01] A. V. Knyazev. Toward the optimal preconditioned eigensolver: locally optimal block pre-
conditioned conjugate gradient method. SIAM J. Sci. Comput., 23(2):517–541, 2001. Copper
Mountain Conference (2000).
[Kuc01] P. Kuchment. The Mathematics of Photonic Crystals, chapter 7, pages 207–272. SIAM,
Philadelphia, PA, 2001.
[Mię10] A. Międlar. Inexact Adaptive Finite Element Methods for Elliptic PDE Eigenvalue Problems.
PhD thesis, TU Berlin, 2010.
[MM11] V. Mehrmann and V. Miedlar. Adaptive computation of smallest eigenvalues of self-adjoint
elliptic partial differential equations. Numer. Linear Algebra Appl., 18(3):387–409, 2011.
[Par98] B. N. Parlett. The symmetric eigenvalue problem. SIAM, Philadelphia, PA, 1998.
[Sim13] S. H. Simon. The Oxford Solid State Basics. Oxford University Press, New York, 2013.
