OPTIMAL DESIGN PROBLEMS IN THIN-FILM AND DIFFRACTIVE OPTICS
By
Yuliang Wang

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
Applied Mathematics - Doctor of Philosophy
2013

ABSTRACT
OPTIMAL DESIGN PROBLEMS IN THIN-FILM AND DIFFRACTIVE
OPTICS
By
Yuliang Wang
Optical components built from thin-ﬁlm layered structures are technologically very important. Applications include but are not limited to energy conversion and conservation,
data transmission and conversion, space technology, imaging and so on. In practice these
structures are deﬁned by various parameters such as the refractive-index proﬁle, the layer
thickness and the period. The problem to ﬁnd the combination of parameters which yield
the spectral response closest to a given target function is referred to as optimal design. This
dissertation considers several topics in the mathematical modeling and optimal design of
these structures through numerical optimization.
A key step in numerical optimization is to deﬁne an objective function the measures the
discrepancy between the target performance and that of the current solution. Our ﬁrst topic
is the impact of the objective function, its metric in particular, on the optimal solution and
its performance. This is done by numerical experiments with diﬀerent types of antireﬂection
coatings using two-material multilayers. The results conﬁrm existing statements and provide
a few new ﬁndings, e.g. some speciﬁc metrics can yield particularly better solutions than
others.
Rugates are optical coatings presenting continuous refractive-index proﬁles. They have
received much attention recently due to technological advances and their potential better
optical performance and environmental properties. The Fourier transform method is a widely
used technique for the design of rugates. However, it is based on approximate expressions

with strict assumptions and has many practical limitations. Our second topic is the optimal
design of rugates through numerical optimization of objective functions with penalty terms.
We found solutions with similar performance and novel solutions by using diﬀerent metrics
in the penalty term.
Existing methods used only local basis functions such as piece-wise constant or linear
functions for the discretization of the refractive-index proﬁle. Our third topic is the use global
basis functions such as sinusoidal functions in the discretization. A simple transformation is
used to overcome the diﬃculty of bound constraints and the result is very promising. Both
multilayer and rugate coatings can be obtained using this method.
Diﬀraction gratings are thin-ﬁlm structures whose optical properties vary periodically
along one or two directions. Our ﬁnal topic is the optimal design of such structures in
the broadband case. The objective functions and their gradient are obtained by solving
variational problems and their adjoints with ﬁnite element method. Interesting phenomena
are observed in the numerical experiments. Limitations and future work in this direction are
pointed out.

ACKNOWLEDGMENTS

I would like to express my deep gratitude and appreciation to my dissertation advisor, Prof.
Gang Bao for his support, understanding, encouragement, and patient guidance over the
long term of my graduate study.
I would also like to thank my committee members Prof. Jianliang Qian, Prof. Chang Yi
Wang and Prof. Zhengfang Zhou for their valuable kindness and help.
It is a great pleasure to thank my colleagues and friends for their solicitude and helpful
suggestions.

iv

TABLE OF CONTENTS

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

8
8
10
13

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

19
21
24
29

Chapter 4 Rugate Coatings . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Smooth by Penalty Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
36

Chapter 5 Global Basis Functions . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40
40
42

Chapter 6 Diﬀraction Gratings .
6.1 Model and Forward Solver . .
6.2 Adjoint Problem and Gradient
6.3 Numerical Experiments . . . .

.
.
.
.

45
45
50
52

Chapter 7 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Chapter 2 Basic Theory for Thin-ﬁlm Optics
2.1 Maxwell’s Equations . . . . . . . . . . . . .
2.2 Model Equations for Layered Medium . . . .
2.3 Reﬂectance and Transmittance . . . . . . .
Chapter 3 Norm of Objective Functions
3.1 Formulation and Algorithms . . . . . .
3.2 Numerical Experiments and Discussion
3.3 Conclusion . . . . . . . . . . . . . . . .

BIBLIOGRAPHY

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

. . . . . . . . . . . . .

v

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

. . . . . . . . . . . . . . . . .

62

LIST OF FIGURES

Figure 2.1

Problem geometry for multilayer coatings. . . . . . . . . . . . . . . .

15

Figure 3.1

(For interpretation of the references to color in this and all other
ﬁgures, the reader is referred to the electronic version of this dissertation.) Dependence of the average and maximal power reﬂectance
on the norm of objective functions for three types of AR coatings.
Top to bottom: broadband, omnidirectional and broadband-andomnidirectional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

Refractive-index proﬁles and power of reﬂection in the omnidirectional case for p = 98 (top) and p = 99 (bottom). . . . . . . . . . . .

29

Figure 3.2

Figure 3.3

Refractive-index proﬁles and reﬂected power for the broadband case
at p = 1, 100 and p = ∞ (top to bottom). . . . . . . . . . . . . . . .

30

Refractive-index proﬁles and reﬂected power for the omnidirectional
case at p = 1, 100 and p = ∞ (top to bottom). . . . . . . . . . . . .

30

Refractive-index proﬁles and reﬂected power for the broadband-andomnidirectional case at p = 1, 100 and p = ∞ (top to bottom). . . .

31

Figure 4.1

Problem geometry for rugate coatings. . . . . . . . . . . . . . . . . .

34

Figure 4.2

Refractive-index proﬁles (left) and reﬂectance curves (right) of optimal solutions for a ramp ﬁlter. Top to bottom: α = 0, 0.01, 0.1. . . .

37

Refractive-index proﬁle (left) and reﬂectance curve (right) of the optimal solution for a ramp ﬁlter with α = 0.01, p = 2, q = 10. . . . . .

38

Figure 4.4

Dependence of the solution on N . Top to bottom: N = 64, 128, 256.

38

Figure 5.1

Snapshots of the refractive-index proﬁle (left) and reﬂectance (right)
at 10, 50, and 200 iterations (top to bottom) during the optimization
of a broadband AR coating parameterized by sinusoidal functions. .

43

Problem geometry for 1D diﬀraction gratings. . . . . . . . . . . . . .

46

Figure 3.4

Figure 3.5

Figure 4.3

Figure 6.1

vi

Figure 6.2

Figure 6.3

Figure 6.4

Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection coating with the homogeneous proﬁle as the initial guess.
Top: initial guess; Bottom: optimized. . . . . . . . . . . . . . . . . .

54

Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection coating with the rectangular proﬁle as the initial guess.
Top: initial guess; Bottom: optimized. . . . . . . . . . . . . . . . . .

56

Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection coating with the triangular proﬁle as the initial guess.
Top: initial guess; Bottom: optimized. . . . . . . . . . . . . . . . . .

58

vii

Chapter 1
Introduction
For almost any practical applications of optics we will encounter a situation when light
travels from one medium to another. In the optical community, one calls the input light the
incident light. The medium from which the incident light propagates is called the incident
medium and the other medium is called the substrate. The behavior of an incident light
is completely determined once the properties of the incident medium and the substrate are
given. An optical coating is deﬁned as an additional layer of medium inserted between the
two media in order to alter the behavior of incident light and achieve certain optical eﬀects,
such as reduction or enhancement of reﬂections. This work is concerned with the optimal
design of optical coatings by numerical methods. The coatings considered are divided into
two major types: classical thin-ﬁlm coatings and diﬀraction gratings.
A thin-ﬁlm coating consists of stratiﬁed media so that its optical properties vary in the
direction normal to the interface but are invariant in the plane parallel to the interface.
Probably the easiest way to obtain a thin-ﬁlm optical coating is by putting a drop of oil on
the surface of water. The interference eﬀects at the interfaces between air, the oil ﬁlm and
water causes variation of reﬂectance with respect to wavelength and forms the colors as we
see. The non-uniformity of the color is due to the non-uniformity of the thickness of the
ﬁlm and the fact that interference eﬀects changes with the thickness. A similar example is
soap bubble, in which case the coating medium is soap water and air plays the role of both

1

incident medium and substrate.
The coatings just mentioned have no practical application except producing beautiful
color patterns. For a coating to be useful, one need have control over its parameters. The
simplest possible coating would be consisting a single homogeneous layer. From the basic
optical theory for stratiﬁed media, such a coating will reduce the reﬂectance to zero if its
refractive index is

√

na ns and its optical thickness is one fourth of the wavelength. Here

na , ns denotes the refractive indices of the incident medium and substrate respectively.
The optical eﬀects which can be achieved by a single homogeneous layer is very simple and
limited. Mathematically speaking, this is because such a coating has only two parameters
that can be tuned: its refractive index and thickness1 . To achieve more complex eﬀects,
one need to increase the number of parameters of the coating. The most intuitive way to
accomplish this is to use a stack of homogeneous layers. The coatings obtained this way are
called multilayer coatings, often abbreviated “multilayers”.
Multilayers have been successfully designed and produced to fulﬁll a wide range of purposes. Today multilayer coatings are employed in most antireﬂection coatings, dielectric
mirrors, bandpass ﬁlters and various other optical devices. These devices have wide applications in many scientiﬁc and industrial areas such as astronomy, biology, medicine, energy
conversion and conservation, data conversion and transmission, space technology, imaging,
photography and so on.
The determination of the optical properties of a given multilayer is regarded as a forward
problem. The optical properties just mentioned can refer to the relationship of the amplitude,
1 The

number of parameters will be three if permittivity and permeability can be tuned
independently. This may be achievable with modern technology, but it is out of the scope
of this work. Moreover we assume the real and imaginary parts of the refractive index can
not be changed independently.

2

phase, angle, energy and polarization states of the reﬂected, transmitted and absorbed light
on those of the incident light. The forward problem is well-posed in the sense that the
solution is existent, unique and stable with respect to the input parameters. The theory and
techniques for solving the this forward problem is relatively complete. However, the accuracy
and eﬃciency of the solution methods are still very important issues since the optimal design
problems usually requires solving a large number forward problems. In chapter we shall
summarize the basic theory and techniques pertaining to the optimal design problems.
Conversely, the determination of multilayers that yield desired optical properties is regarded as an inverse problem. Unlike forward problems, inverse problems are usually illposed, namely the solutions may not be existent, unique or stable with respect to the input
parameters. In particular, inverse problems for practical optical coatings usually have no
exact solution. Consequently we have the take the second best and look for approximate
solutions whose output are as close as possible to the desired ones. This is termed as optimal
design problems.
Optimal design problems can be considered as generalized inverse problems and share
some common properties with the latter, such as the non-uniqueness of solutions and their
high sensitivity with respect to the target performance. In general these properties are
considered as disadvantages for an inverse problem since its objective is to determine the
solution, and the solution should be as accurate as possible. However, they could be advantages for an optimal design problem since its objective is to ﬁnd solutions yielding the
best performance. Non-uniqueness means we can have diﬀerent designs giving the same
performance, which is beneﬁcial since we have an opportunity to select the best feasible one.
High sensitivity of the solution with respect to the performance implies low sensitivity of the
performance with respect to the solution, an advantage in practice indeed.
3

Before computers became powerful and cheaply available, multilayers are designed by
analytic techniques such as the vector method, the Herpin index, the method of eﬀective
interfaces, the Smith’s chart and so on. A complete account of these and other similar
methods can be found in the books by Furman and Tikhonravov [1], Macleod [2] and the
references cited therein. The analytic methods are intrinsically limited to small number of
layers and simple target speciﬁcations. In addition, they can not be considered as optimal
design methods. However, they are still useful today since they can provide insight and
intuitive explanations which are generally lacking in numerical methods. In particular, the
solutions obtained by analytic methods can be used as initial guess for numerical optimization
methods.
The advances of modern computers boosts the development of numerical methods, numerical optimization methods in particular. Baumeister [3] is probably the ﬁrst one to
employ numerical optimization methods in the design of optical coatings. Since then many
standard algorithms have been investigated for their applications in the optimal design of
multilayers. For example, Dobrowolsk [4] applied and compared ten algorithms, including
both deterministic and stochastic types. There exist various algorithms developed speciﬁcally for multilayers, too. One of the most notable is the so called needle method, originally
proposed by Tikhonravov [5] and then developed and matured by himself and others. This
method has been shown eﬀective in the design of multilayers for a wide range of applications.
However, none of the methods can claim to be the best for all cases due to the variety of the
applications. In this work we shall not consider the problem of ﬁnding or devising the best
algorithm for all design problems.
A key step in solving a practical problem by numerical optimization is to deﬁne an
objective function to be minimized. The objective function should measure the discrepancy
4

between the target performance and that of the current solution. The metric in which the
discrepancy is measured will certainly aﬀect the solutions and their performance. Chapter is
devoted to our study on this subject. We will demonstrate the impact of the metrics through
numerical examples and observe both anticipated and unanticipated phenomena. Practical
implication and value of the observations will be pointed out.
As mentioned earlier, using a stack of homogeneous layers is one way to increase the
degree of freedom of a thin-ﬁlm structure so as to achieve complex and powerful optical
eﬀects. Another approach to accomplish the goal is to use a single inhomogeneous, mostly
continuous layer. The coatings obtained this ways are called rugate coatings, often abbreviated “rugates”. Traditionally rugates refer to coatings whose refractive-index proﬁles are
constructed from sinusoidal and other elementary functions and which are much thicker than
their multilayer counterparts. Oftentimes this name is used interchangeably with “gradientindex coatings” and “inhomogeneous coatings”, both of which refer to coatings with arbitrary
continuous refractive-index proﬁles and similar thickness compared to multilayers. We shall
pose no restrictions on the thickness of coatings and use “rugates” for all coatings with
continuous refractive-index proﬁles.
Rugates have received much attention recently due to their better environmental properties, less sensitivity to production errors and potential better optical performances in complicated situations. The Fourier transform method is a widely used and powerful technique for
the design of rugates. However, it is based on approximate expressions with strict assumptions and is restricted to single incident angle and polarization state. Another approach is to
divide the coating into a large number of thin homogeneous layers and optimize an objective
function with respect to the refractive indices of all layers. Because of the maximum principle of thin-ﬁlm optics [6], the coatings obtained this ways generally resemble multilayers
5

or possess many discontinuous points if the usual objective functions are applied. In order
to force smoothness to the solutions, we add penalty terms to the usual objective functions.
In chapter we study these methods and obtain similar solutions as other methods and new
ones never seen before by using diﬀerent norms in the penalty terms.
Existing optimization methods for the optimal design of thin-ﬁlm coatings have used only
local basis functions (step and hat functions) for the parameterization of the refractive-index
proﬁle. Global basis functions such as sinusoidal functions have been used only in heuristic
design methods. In chapter we propose a method to use a global basis functions in the
design of optical coatings. The problem of using global basis functions in this context is
that there is no direct way to translate the bound constraints on the refractive-index proﬁle
to its spectral coeﬃcients. We overcome this diﬃculty by a simple transformation of the
bounded functions to unbounded ones. This way the truncated series expansion can be used
to approximate the coating proﬁle and unconstrained optimization algorithms can be used
on the series coeﬃcients. An interesting feature of this method is that we can obtain both
rugates and multilayers with very similar performances.
So far we considered only optical coatings built from layer structures. The refractive
index of such structures vary along one direction and is homogeneous over the other two.
By allowing the refractive index of the structure to vary periodically along one or both of
the other two directions one obtains singly or doubly periodic structures, often called 2D
or 3D diﬀraction gratings, respectively. Dobson studied the optimal design of wide-angle
antireﬂective diﬀraction gratings [7] and blazed diﬀraction gratings [8]. Following the spirit
of these work, we study the optimal design of broadband antireﬂective diﬀraction gratings
in chapter . The value and gradient of the objective function was obtained by solving a
2D PDE with nonlocal transparent boundary conditions [9] and the corresponding adjoint
6

equation [7], respectively.

7

Chapter 2
Basic Theory for Thin-ﬁlm Optics
The basic theory for thin-ﬁlm optics can be found from a few texts, such as Born and
Wolf [10], Furman [1] and Macleod [2]. In this chapter we summarize the parts of the theory
that are essential and necessary for our work. We start from the Maxwell’s equations and
derive the formulas to compute the reﬂectance, transmittance and their derivatives with
respect to the layer thicknesses and refractive indices of a multilayer system.

2.1

Maxwell’s Equations

The behavior of light, and electromagnetic ﬁeld in general, is described by the Maxwell’s
equations (in the Gaussian system of units):

4π
1 ∂D
=
j,
c ∂t
c
1 ∂B
×E +
= 0,
c ∂t

×H−

(2.1a)
(2.1b)

· D = 4πρ,
· B = 0,
where

8

(2.1c)
(2.1d)

E

electric vector

H

magnetic vector

D

electric displacement

B

magnetic induction

j

electric current density

ρ

electric charge density

c

speed of light in vacuum

In the ﬁeld of thin-ﬁlm coatings the source terms j and ρ are assumed to be zero. In addition,
the other four quantities E, H, D, B must be supplemented by material equations for a unique
determination of the ﬁeld vectors. They are rather complicated in general. In this work we
assume the body is moving at slow speed, the material is isotropic and linear, and the ﬁeld
is time harmonic. Under these assumptions the material equations are give by

D = εE,
B = µH.
where ε and µ are called the relative permittivity and permeability respectively.
We assume the ﬁeld is time harmonic and write

E = E exp(−iωt),
H = H exp(−iωt),
where ω is the angular frequency of the ﬁeld.

9

Applying all the above assumptions to (2.1) we arrived at these equations

× E − ikµH = 0,

(2.2a)

× H + ikεE = 0,

(2.2b)

where k = ω/c is the wavenumber in vacuum. Notice k = 2π/λ if λ is the wavelength in
vacuum.
Finally we have the following boundary conditions at a surface of discontinuity:

n12 × (E2 − E1 ) = 0,

(2.3a)

n12 × (H2 − H1 ) = 0,

(2.3b)

where n12 is the normal vector of the surface from region 1 to region 2 and Ej , Hj are the
ﬁeld vectors in region j respectively.

2.2

Model Equations for Layered Medium

For the sake of simplicity we also assume the material is nonmagnetic, namely the permeability µ = 1 everywhere, so that the problem geometry is characterized completely by the
permittivity ε. A thin-ﬁlm coating consists a layered medium of ﬁnite depth sandwiched

10

between two homogeneous media of inﬁnite depth. In other words, we have



 ,
a
z ≤ 0,





ε(x, y, z) =
 (z), 0 < z < L,







 s,
z ≥ L.
The medium occupying the region z ≤ 0 is called the incident medium and the medium
occupying the region z ≥ L is called the substrate. The incident medium is where the source
light (incident light) comes from. We assume the incident light is a plane wave, namely
the ﬁeld vectors E, H have spatial dependence exp(ik · r) where k is a wave vector with
√
|k| = k εa and r = (x, y, z) is the spatial vector. Notice k is the direction of propagation of

the light. Due to the invariance of ε along the x, y axes, any plane wave solution of (2.2) can
be decomposed into two components called transverse-electric (TE) and transverse magnetic
(TM) modes, respectively.
The plane formed by the propagation direction k and the z-axis is termed the incident
plane. Since ε is invariant along both x and y axes, we assume the incident plane is the yz
plane without loss of generality. In the case of TE mode, the electric ﬁeld is perpendicular
to the incident plane. So E has only the x-component diﬀerent from zero. Let us write
E = (u, 0, 0) and substitute it into (2.2). Expanding the

× operators we ﬁnd H lies in

the yz plane and is independent on x. Hence u is also independent on x and we can write
u = u(y, z). The resultant equation for u is the following:

(∆ + k 2 ε)u = 0,

11

(2.4)

where ∆ = (∂yy + ∂zz ). By the variable separation method we ﬁnd

u(y, z) = u(z) exp(ikαy),

where α is a constant to be determined. Substituting this to (2.2) yields

H = (0, −v(z), αu(z)) exp(ikαy),
where v and u are solutions to the following system of linear ordinary diﬀerential equations:

du
= ikv,
dz
dv
= ik(ε − α2 )u.
dz

(2.5a)
(2.5b)

To determine the value of α, consider the wave in the incident medium. If the incident angle
is θI , then the wavevector
√
k = k εa (0, sin θI , − cos θI )
and
√
u(y, z) = exp(ik εa )(y sin θI − z cos θI ).
Hence α =

√
εa sin θI .

In the case the TM mode, the H ﬁeld is perpendicular to the incident plane and we can
write H = (u(z), 0, 0). Similar to the case of TE mode, we ﬁnd

α
E = 0, v(z), − u(z) exp(ikαy),
ε
12

and u, v are solutions to the following system of linear ordinary diﬀerential equations

du
= ikεv,
dz
dv
α2
= ik 1 −
dz
ε

2.3

(2.6a)
u.

(2.6b)

Reﬂectance and Transmittance

Now we derive formulas for the reﬂectance and transmittance for a layered structure. Let us
assume TE mode for the time being. For the sake of convenience let us normalize the ﬁeld
so that the electric ﬁeld at the substrate interface u(L) = 1. Denote uI , vI , uR , vR to be the
ﬁeld of the incident wave and reﬂected wave respectively. By (2.5) we obtain

vI (0) = qa uI (0),

(2.7)

vR (0) = −qa uR (0),
where qa =

εa − α2 . By the jump conditions (2.3) we obtain the continuity of u, v at the

interface z = 0. As a result, we have

u(0) = uI (0) + uR (0),
v(0) = vI (0) + vR (0) = qa (uI (0) − uR (0)).
Solving the above system for uI (0) and uR (0) we obtain the following
qa u(0) + v(0)
,
2qa
qa u(0) − v(0)
uR (0) =
,
2qa
uI (0) =

13

The amplitude reﬂectance r is deﬁned to be the ratio of the tangential component of the
electric ﬁeld of the reﬂected wave and incident wave. Hence we obtain r in terms of the total
ﬁeld at the incident interface z = 0:
qa u(0) − v(0)
u (0)
=
.
r= R
uI (0)
qa u(0) + v(0)

(2.8)

Similarly we obtain the amplitude transmittance t in terms of the total ﬁeld at the
substrate interface z = L:
t=

u(L)
2qa
=
.
uI (0)
qa u(0) + v(0)

(2.9)

Now consider the special case when the coating layer is homogeneous, namely we assume
ε(z) = ε, 0 < z < L. For the sake of convenience we introduce the refractive index n =

√
ε.

Eliminating v from (2.5) yield the Helmholtz equation for u:
d2 u
+ k 2 (n2 − α2 )u = 0.
dz 2

(2.10)

Denote β 2 = k 2 (n2 − α2 ), then the general solution of the above equation is the following
u(z) = c1 exp(iβz) + c2 exp(−iβz).

From the ﬁrst equation in (2.5) we obtain

v(z) = q(c1 exp(iβz) − c2 exp(−iβz)),
where we denote q =

β
. Substituting z = 0 and z = L in the above solutions respectively and
k

eliminating c1 , c2 , we obtain the following equation relating the ﬁeld values at the substrate
14

x
Incident medium
n = na

z
ni

···

θ

···

ti

Incident light

y
Substrate
n = ns

Figure 2.1: Problem geometry for multilayer coatings.
interface to those at the incident interface:




i
cos(βL)
− sin(βL)   u(L) 
 u(0)  
q


=



v(0)
v(L)
−iq sin(βL)
cos(βL)




Let u(L) = 1, then v(L) = qs =

(2.11)

εs − α2 from continuity. Plugging these into (2.11)

we obtain u(0), v(0). Together with (2.8) and (2.9) we can write the expressions for the
amplitude reﬂectance r and amplitude transmittance t.
Now consider a multilayer coating with layer thicknesses {t1 , · · · , tn } and refractive indices {n1 , · · · , nN } (See Figure 2.1). Combining (2.11) for all layers we obtain




N
 u(0)  
 u(L) 
Mj  

=
,
j=1
v(0)
v(L)

where L =





n
j=1 tj ,

βj =


i
− sin(βj tj ) 
 cos(βj tj )
qj
,
Mj = 


−iqj sin(βj tj )
cos(βj tj )


(2.12)

βj
n2 − α2 and qj =
. Then we can calculate r and t as in
j
k

the single layer case. This technique for the computation of reﬂectance and transmittance
of a multilayer is called the matrix method and will be used through out this work for the

15

computation of the objective functions.
In order to use gradient based optimization algorithms such as quasi-newton methods, one
need the gradient of the objective functions, which in turn are computed from the gradient
of r and t. From (2.8) and (2.9) we see the gradient of r and t can be expressed in terms of
the gradient of u(0) and v(0), which can be easily obtained from (2.12) since each parameter
tj0 and nj0 appears only in Mj0 . Therefore we have






j0 −1



dMj0
d  u(0)  
Mj 

=
dtj0
dtj
j=1
v(0)

 

j0 −1
dMj0
d  u(0)  
Mj 

=
dnj0
dnj
j=1
v(0)



 u(L) 

Mj  
,
j=j0 +1
v(L)



N
 u(L) 

Mj  
.
j=j0 +1
v(L)


N

(2.13a)

(2.13b)

The computational cost for r and t is O(N ) due to multiplication of N 2 × 2 matrices.
The computational cost for each partial derivatives is also O(N ) and the computational
cost for the gradient will be O(N 2 ) if done straightforwardly using (2.13). The cost can be
reduced to O(N ) once we notice the matrix products
j0 −1
j=1

N

Mi ,

Mj .
j=j0 +1

can be computed progressively, once in forward direction and another in the backward direction, as noted by, e.g. Verly [11]. In addition, we utilized the Armadillo C++ linear algebra
library [12] for computation of matrix multiplications due to its exceptional eﬃciency in
handling small matrices.
Finally, in many practical situations we are concerned with only the energy carried in
the reﬂected and transmitted waves. The magnitude and direction of the energy carried by
16

a time-harmonic EM wave is given by the time-averaged Poynting vector

S=

c
(EH).
8π

Using out notations the energy carried by the wave along the z-axis is given by

Sz =

c
(uv).
8π

From (2.7) we ﬁnd the energy carried by the incident and reﬂected wave along the z-axis is
given by
c
(qa )|uI (0)|2 ,
8π

−

c
(qa )|uR (0)|2 ,
8π

respectively. The power reﬂectance R is deﬁned by the ratio of the energy carried by the
reﬂected wave along the z-axis to that of the incident wave. Thus we have

R = |r|2 .

(2.14)

Similarly we can derive the formula for the power transmittance:

T =

(qs ) 2
|t| .
(qa )

(2.15)

The gradient of R and T with respect to layer thicknesses and refractive indices are computed

17

easily from those of r and t:

R = 2 (r r),
T =2

(qs )
(t t).
(qa )

18

(2.16a)
(2.16b)

Chapter 3
Norm of Objective Functions
In this chapter we consider the optimal design of two-material multilayer AR coatings. The
most common approach to this problem is numerical optimization. Baumeister [3] seems
to be the ﬁrst one to employ this approach. Since then people have investigated various
existing optimization algorithms and developed many new ones. For example, ten diﬀerent
algorithms of deterministic and stochastic types were tested and compared in [4]. The needle
method proposed in [5] is one of the most powerful algorithms. It has been successfully used
in the design of not only antireﬂection but also other complicated optical coatings.
Besides the conventional stratiﬁed thin-ﬁlm coatings, numerical methods have been studied and successfully applied to the optimal design of other thin-ﬁlm structures, such as
diﬀraction gratings with speciﬁed diﬀraction patterns [13, 14, 15, 9, 8] and with omnidirectional antireﬂection properties [7], guided-mode grating resonant structures [16, 17, 18] and
nonlinear diﬀraction gratings [19].
A key step in numerical optimization is to deﬁne a merit function that measures the
discrepancy between the target and and solution. The metric in which the discrepancy is
measured will certainly aﬀect the optimal solution and its performance. The most commonly
used metric is of the least-square type. Others include the sum of residuals or the max of the
residuals. For broadband AR coatings, the merit functions corresponding to the previously

19

mentioned metrics are given by

1
M
1
2
[R(ki )]2
,
M

i=1

1
M

M

R(ki ),
i=1

max R(ki ),

1≤i≤M

where R(ki ) denotes the reﬂectance at wavenumber ki and {ki : 1 ≤ i ≤ M } is a set of
uniformly distributed sampling points in the spectral region of interest.
Mathematically speaking, these can be considered as discrete approximations of the
L2 , L1 and L∞ norms of the function R(k) over the spectral region. Naturally we can
also consider using Lp norms for other positive integers p. It was mentioned in [20] that
greater p will force solutions which tend to equalize the residuals. However, no theoretical
or numerical substantiations have been made.
In this chapter we conduct quantitative studies on the impact of p on the optimal design
of antireﬂection coatings . This is accomplished by comparing the average and maximal
values of the reﬂected power over the spectral and/or angular regions of interest through
numerical experiments. These speciﬁc values are chosen because they are the most important
for many practical applications. The average value, up to a weight function, measures the
average power reﬂected by the interface. Thus it is more important than other quantities in
situations when maximal power absorption is mostly desired, e.g. solar cell panels. On the
other hand, imaging devices such as glasses and camera lenses should be designed so that
the reﬂected power is more evenly spread and the maximal value is as small as possible.

20

3.1

Formulation and Algorithms

Without loss of generality, the AR coating will be assumed to be of broadband type in this
section. Formulations for omnidirectional and broadband-and-omnidirectional types are similar. The spectral region of interest will be denoted by an interval [kL , kH ] of wavenumber or
an interval [λL , λH ] of wavelength. Wavenumber is better suited in mathematical formulation while wavelength seems to be more natural in the description of practical applications.
Let {ki : 1 ≤ i ≤ M } be a set of uniformly distributed points in [kL , kH ] with k0 = kL and
kM = kH and deﬁne a class of merit functions {fp : 1 ≤ p ≤ ∞} by

1



p
 1 N



, if p < ∞,
[R(ki )]p

fp =  N i=1





 max R(k ),

if p = ∞,
i

(3.1)

1≤i≤N

where R is reﬂectance as a function of wavenumber. Dispersion of all materials are
ignored in the calculation of fp and its gradient, although it can be incorporated into our
model through simple transformations as in [21].
The value of fp can be computed directly from the values of R(ki ) for all p. The partial
derivatives of fp can be computed directly from the values of R(ki ) and their partial derivatives for p < ∞. The function f∞ , on the other hand, is not diﬀerentiable. In order to use
gradient-based optimization algorithms, we transform the original optimization problem

min max R(L, ki )
L 1≤i≤N

21

into an equivalent form:

min s subject to R(L, ki ) − s ≤ 0 ∀i,
L

(3.2)

where L = {t1 , · · · , tN } is the design variable. The equivalence of these two problems is
shown in e.g. [22] and it is the approach adopted in [23]. Note that the objective function
and all the constraint functions in (3.2) are diﬀerentiable. Diﬀerent approaches have been
proposed in [24, 25]. They are more complicated but can be used to design other types of
coatings.
Now that the optimization problems are completely setup, we need to choose algorithms.
Although the needle method could be the best candidate, it is rather complicated compared
with straightforward optimization algorithms. Since our objective is not to ﬁnd the best
algorithm but to study the impact of the metrics, we choose to utilize existing algorithms in
the open source library NLopt [26]. All the local gradient-based algorithms in this library
can handle cases for p < ∞. Among them we ﬁnd the low-storage BFGS algorithm described
in [27] to be particularly eﬃcient for our test problems. If it fails to converge, we can always
turn to more robust but less eﬃcient ones. The only algorithm in this library that can
handle nonlinear constraints is a sequential quadratic programming method based on the
implementation in [28], so it is the one we shall use for the case p = ∞.
For two-material coatings, the refractive indices ni take a low and high value alternatively.
They are ﬁxed up to the value of the ﬁrst layer. Thus the design variables on which we do
optimization are layer thicknesses ti . For both practical and computational considerations,
we need deﬁne lower and upper bounds for the design variables. First of all, we need set a
lower bound, which must be non-negative in order to avoid nonphysical solutions. In fact, it

22

will be set to 0 in our numerical experiments. However, zero-thickness and thin layers could
occur in the solution. Due to limited accuracy of equipment and sensitivity to manufacturing
errors, thin layers should be avoided in the design. This can be achieved by specifying a
positive value for the lower bound in the beginning or by a process named “consolidation” [29,
11], which eliminates thin layers during or after the optimization. There exist two possible
approaches for this process: erase thin layers or merge them with neighboring thicker layers.
Existing works did not mention which method is better and which one was used in their
works. We tested both approaches and found the latter usually yields better results than
the former. Our algorithm for the consolidation process is the following:
set thresholds Ti for 1 ≤ i ≤ N ;
for i = 2 to M do
if ti < Ti then
ni ← ni−1 ;
end
if t1 < T1 then
n1 = n2 ;
end
end
combine neighboring layers with the same refractive index
Algorithm 1: Consolidation of multilayer coatings.

The procedure “optimization → consolidation” is iterated until no thin layers remain.
The thresholds ti are user input values that specify the least thicknesses to be present in the
coating and can be chosen independently for each layer. Obviously they should not be too
large, since otherwise the obtained solution will be far from optimal. They can be very small
so that only the layers with zero thicknesses are consolidated. Interestingly, by using suitable
values we obtain coatings with much less layers but similar or even better performances.
An upper bound is not needed unless it is a required feature of the coating since our
algorithms always generate solutions without too thick layers.
23

3.2

Numerical Experiments and Discussion

Besides the merit functions, the lowest achievable residual reﬂectance of an antireﬂection
coating depends on many other parameters such as
• refractive indices of the coating material and the ambient media
• total number of layers
• total optical thickness
• spectral and/or angular regions
• state of polarization
An explicit formula for this dependence is not available, but some empirical expressions
have been proposed in [30, 31, 32, 33]. For the purpose of our numerical experiments, these
parameters will be held ﬁxed when p is varied.
In all the following experiments, the refractive indices of the incident and substrate media
are 1.0 and 1.52, corresponding to air and glass, while the low and high refractive indices of
the two materials used in the coating are 1.45 and 2.35, respectively. The number of layers is
set to 20 in the beginning and will be reduced in the ﬁnal solutions due to the consolidation
process.
Initial guesses play an important role for local optimization methods. We use a quarterwave stack at a ﬁxed wavelength λc . That is, we set ti = λc /(4ni ) for each i. In the case
of broadband AR coatings, good results can be obtained if λc is slightly less than λL as
suggested in [34]. According to our numerical experiments, smaller λc will cause more layers
to have zero thicknesses during the optimization. This leads to less number of layers and

24

worse performance in the ﬁnal solution. If λc is too small, all layers with refractive index nH
will be eliminated and the ﬁnal solution consists of a single layer with refractive index nL .
On the other hand, the solution will be far from optimal if λc > λL . The best choice should
be made by a trial-and-error method and is dependent on the speciﬁc problem. However, our
numerical experience shows that the choice λc = 0.8λL can produce good results consistently
for many diﬀerent combinations of parameters.
The above observations are true not only for the broadband type but also for the omnidirectional (with λ0 playing the role of λL ) and broadband-and-omnidirectional types. In
all the results shown in the sequel we used the initial guesses with λc = 0.8λL .
As mentioned in section 3.1, the thickness thresholds ti controls a balance between the
performance and the minimal layer thickness of the ﬁnal coating. We used ti = λL /(60ni )
for all the experiments.
We now present results for the following three AR coatings
(a) broadband: normal incidence over spectral region [λL , λU ] = [400, 1200]
(b) omnidirectional: at ﬁxed wavelength λ0 = 100 over angular region [θL , θH ] = [0◦ , 90◦ ]
(c) broadband-and-omnidirectional: over the spectral region [λL , λU ] = [400, 1200] and angular region [θL , θH ] = [0◦ , 90◦ ].
The merit functions for the three cases are deﬁned as the following respectively:

1



p
 1 M


p

[R(ki )]
, if p < ∞,

fp =  M i=1





 max R(k )I(k ),

if p = ∞,
i
i
1≤i≤M

25


1



p
 1 N


p

[R(θi ) cos(θi )]
, if p < ∞,

fp =  N i=1





 max R(θ ) cos(θ ),

if p = ∞,
i
i
1≤i≤N


1



M N
p


1


[R(ki , θj ) cos(θj )]p
, if p < ∞,

fp =  M N i=1 j=1







max
R(ki , θj ) cos(θj ),
if p = ∞,
1≤i≤M,1≤j≤N

where {k1 , · · · , kM } and {θ1 , · · · , θN } are uniformly distributed sampling points in [kL , kU ]
and [θL , θU ], respectively. Notice our deﬁnitions of merit functions for the omnidirectional
and broadband-and-omnidirectional cases are diﬀerent from conventional ones by a factor
of cos(θi ). This factor is included so that the merit functions represent the reﬂection in
the sense of total power. This diﬀerence is not essential for our discussion since it can be
regarded as a weight function. The observations and conclusions in the sequel are still valid
if this term is dropped out.
Notice that the average and maximal power of reﬂection Pave and Pmax are approximated
by f1 and f∞ respectively. In fact, they are computed by the same formula except with with
a larger number of sampling points in our experiments. Intuitively the values of Pave , Pmax
should be the lowest when f1 , f∞ are used as the merit functions and should increase and
decrease as p increases when fp are used as merit functions. These are generally the case as
seen from Figure 3.1, where these values are plotted against p for 1 ≤ p ≤ 100.
There are several other interesting phenomena associated with Pave and Pmax in Figure 3.1:
• The increase and decrease of the two values slows down quickly as p increases.
26

0.030

Pave
Pmax

0.025
0.020
0.015
0.010
0.005
0
0.016
0.014
0.012
0.010
0.008
0.006
0.004
0.002
0.000
0

20

40

p

60

80

100

80

100

80

100

Pave
Pmax

20

40

p

60

0.14

Pave
Pmax

0.12
0.10
0.08
0.06
0.04
0.02
0

20

40

p

60

Figure 3.1: (For interpretation of the references to color in this and all other ﬁgures, the
reader is referred to the electronic version of this dissertation.) Dependence of the average
and maximal power reﬂectance on the norm of objective functions for three types of AR
coatings. Top to bottom: broadband, omnidirectional and broadband-and-omnidirectional.
27

• The two values becomes rather close at large p in the broadband case but not the other
two. This implies the power of reﬂection will be more evenly spread over the region of
interest in the broadband case than the other two (compare Figures 3.3-3.5).
• The best achievable average and maximal power of reﬂection are much higher ( 10
times) in the broadband-and-omnidirectional case than the other two.
• There exist speciﬁc values of p at which both values experience abrupt changes in the
omnidirectional case. However, this phenomenon is less evident in the broadband and
broadband-and-omnidirectional cases. We believe it is because the merit functions are
smoother in these cases but do not have a theoretical or numerical proof.
The last observation worthies more investigation since it is not as expected as the others.
First we would like to point out the abrupt changes are not caused by the consolidation
process. Our other experiments show they exist even if consolidations are not performed,
although at possibly diﬀerent values of p. At ﬁrst glance it is against the intuition since the
dependence of fp on p is smooth. A small change in p can cause only a small change in fp .
However, fp usually possess a large number of minima and our algorithms are local, so a
small change in fp could drive the solution to a diﬀerent local minimum. For example, the
two values experience abrupt changes from p = 98 to p = 99 and back to their normal values
at p = 100. Figure 3.2 shows the refractive-index proﬁles and distribution of reﬂected power
at p = 98, 99.
Comparing these with the result at p = 100 in Figure 3.4 reveals that the solutions
at p = 98, 100 are the same local minimum but the solution at p = 99 is a diﬀerent one.
The existence of these diﬀerent solutions at speciﬁc values of p is a good thing for coating
designers. This is because one solution may be less sensitive to manufacture errors, possesses
28

1.0
2.35

0.8
0.6
1.45
0.4
2.35
0.2
1.45
0.0
0.0
0

Refractive Index

Reﬂected Power

1.0
0.8 Pave = 0.004, Pmax = 0.007 0.01
0.6
0.00
0.4
0.2 Pave = 0.003, Pmax = 0.006 0.01
0.0
0.00
0.8
500.2 100 0.4 200 250 0.8
150 0.6
300 1.00.0 10 0.2 300.4 50 60 70 80 1.0
0
20
40 0.6
Incident Angle (Degree)
Thickness

Figure 3.2: Refractive-index proﬁles and power of reﬂection in the omnidirectional case for
p = 98 (top) and p = 99 (bottom).
better mechanical properties, or be easier to produce than another depending on the available
deposition technology [35]. For example, the solution at p = 99 can be considered better
than those at p = 98 and p = 100 due to its better optical performance and absence of thin
layers.
Figures 3.3-3.5 show the refractive-index proﬁles and power of reﬂection at p = 1, 100 and
∞. Clearly the power distribution becomes more evenly spread as p increases in all three
cases, conﬁrming the statements in [20]. The solutions for p = 100 and p = ∞ could be very
close, as in the broadband and broadband-and-omnidirectional cases, or rather diﬀerent, as
in the omnidirectional case, where the solution at p = ∞ contains much less layers and has
a small decrease in performance than the solution at p = 100. This can be explained by the
existence of neighboring local minima and the the use of diﬀerent algorithms. Hence it is
worthwhile to study the p = ∞ case since it may lead to diﬀerent, possibly better solutions
than ﬁnite values.

3.3

Conclusion

We have studied the problem of designing AR coatings by numerical optimization methods
with diﬀerent metrics. The merit functions in our approach are discrete approximations of
29

1.0
2.35

Refractive Index

0.8

0.4
1.45
2.35
0.2
0.2
200

0.4
400

0.6
600
Thickness

0.8
800

Reﬂected Power
Pave = 0.006, Pmax = 0.027

0.02
0.01
0.00
0.6 P
Pmax = 0.009 0.02
ave = 0.008,
0.01
0.4
0.00
0.2 Pave = 0.008, Pmax = 0.008 0.02
0.01
0.0
0.00
1.0
0.0
0.2
0.8
1.0
400
600 0.4 800 0.6 1000
0.8

1.45
2.35
0.6

1.45
0.0
0.0
0

1.0

Wavelength

Figure 3.3: Refractive-index proﬁles and reﬂected power for the broadband case at p = 1, 100
and p = ∞ (top to bottom).

1.0
2.35
0.8

1.45
2.35
0.6
0.4
1.45
2.35
0.2

1.45
0.0
0.0
0

Refractive Index

1.0

Reﬂected Power
Pave = 0.001, Pmax = 0.012

0.02
0.01
0.00
0.6
Pave = 0.004, Pmax = 0.007 0.02
0.01
0.4
0.00
0.2 Pave = 0.005, Pmax = 0.011 0.02
0.01
0.0
0.00
0.8
0
20
40 0.6
500.2 100 0.4 200 250 0.8
150 0.6
300 1.00.0 10 0.2 300.4 50 60 70 80 1.0
0.8

Incident Angle (Degree)

Thickness

Figure 3.4: Refractive-index proﬁles and reﬂected power for the omnidirectional case at
p = 1, 100 and p = ∞ (top to bottom).

30

Refractive Index

1.0
2.35
0.8

1.45
2.35
0.6
0.4
1.45
2.35
0.2

Incident Angle

0
20
40
60
80
400
0
20
40
60
80
400

Incident Angle

Incident Angle

1.45
0.0
0.0
0

0
20
40
60
80
400

100

0.2
200

300

0.4
400

500

Thickness

Pave = 0.028,

0.6
600

700

0.8
800

Pmax = 0.137

900

1.0

0.150
0.135
0.120

500

600

700

800

Pave = 0.037,

900

1000

1100

Pmax = 0.073

0.105
0.090
0.075
0.060

500

600

700

800

Pave = 0.031,

900

1000

1100

Pmax = 0.072

0.045
0.030
0.015

500

600

700

800
900
Wavelength

1000

1100

0.000

Figure 3.5: Refractive-index proﬁles and reﬂected power for the broadband-andomnidirectional case at p = 1, 100 and p = ∞ (top to bottom).

31

the Lp norm of the reﬂected power over the region of interest. The impact of p on the optimal
solutions are investigated by comparing the average and maximal power of reﬂection, which
are important factors for practical applications. Based on our numerical experiments with
three types of AR coatings, one should set p to 1 if one wishes to minimize the average value
and to ∞ or a large value (e.g. 100) if one wishes to minimize the maximal value. Moreover,
abrupt changes to the solutions can happen at speciﬁc values of p. This provides coating
designers new opportunities in the search of better solutions.
Finally we would like to point out that the observed eﬀect of metrics on the objective
functions are not limited to multilayers or AR coatings. Similar eﬀects are observed for
rugates and other types of coatings during our numerical investigations.

32

Chapter 4
Rugate Coatings
Rugate coatings have been studied actively by scientists and engineers for more than 50
years. They have received considerable attention recently due to their outstanding optical,
physical and feasibility properties.

4.1

Introduction

Based on the shape of the refractive-index proﬁle, thin-ﬁlm optical coatings are generally
divided into two categories: multilayer coatings and rugate coatings. A multilayer coating
consists of a stack of limited number of homogeneous layers. The refractive-index proﬁle
of such a coating is a step function taking a few number of ﬁxed values. The most common coatings of this type are two-component multilayers, which consist of two materials.
Rugate coatings were initially deﬁned as optical coatings whose refractive-index proﬁles are
constructed from sine and other known functions with explicit expressions. Now they are
usually referred as optical coatings presenting any continuous refractive-index proﬁles. In
this respect the term “rugate” may be used interchangeably with “gradient-index” or “inhomogeneous”. Figure 4.1 shows the geometry for rugate coatings.
Since a continuous function can be approximated by step functions up to any given
accuracy and the practical computation of rugate coatings are commonly conducted through
subdivision of the ﬁlm into a large number of homogeneous layers, the boundaries that
33

x
Incident medium
n = na

z
n = n(z)

θ

y
Substrate
n = ns

Incident light
Figure 4.1: Problem geometry for rugate coatings.

deﬁne multilayer and rugate coatings can become unclear. Those in the intermediate region
can be called quasi-rugate coatings, which include multilayers with intermediate refractive
indices, two-component multilayers with a large number of layers, and coatings possessing
both homogeneous and inhomogeneous layers.
Multilayer coatings have been successfully designed and produced to fulﬁll a wide range
of desired optical properties. But it is still desirable to study rugate coatings because they
are advantageous to multilayers in three aspects [36, 35]. Firstly, rugate coatings can have
better optical properties. In the case of nonabsorbing dielectric coatings at s-polarization
or p-polarization at small incident angles, the maximum principle of thin-ﬁlm optics [6] predicts that the optimal refractive-index proﬁle is a step function taking the lowest and highest
permissible values, i.e. the optimal design is a two-component multilayer. However, the conclusion of this principle becomes unclear or invalid if the coating materials are absorbing,
or if the design target is speciﬁed at wide incident angles and/or both polarization. Under
those circumstances the rugate and quasi-rugate coatings can have superior optical performance comparing with multilayers. Secondly, rugate coatings can have better feasibility
properties. Some multilayer designs are very sensitive to manufacturing errors, especially

34

when there are many thin layers and the total number of layers is large. It has been demonstrated experimentally that many rugate coatings possess better feasibility properties than
their multilayer correspondence. Thirdly, rugate coatings have better physical properties,
including mechanical and thermal stability.
The Fourier transform method [37] is a powerful and widely used method for the design
of rugate coatings. It is based on the approximate expression:
1 ∞
[ln n(x)] eikx dx = Q(k)eiφ(k)
2 −∞
where
• n(x) is the refractive-index proﬁle
• Q(k) is an suitable even function of the target transmittance T (k)
• φ(k) is an odd function to ensure n(x) is real (Note: φ(k) is NOT the phase change of
reﬂection or transmission)
z
• x is twice the optical thickness: x = 2 0 n(s)ds

• k is the wavenumber: k = 2π/λ
The Fourier transform method is very fast when it works. However it has the following
disadvantages and limitations:
• based on approximate equations assuming the following conditions
– low reﬂection
– no absorption (real refractive index)
35

– no dispersion (refractive index is independent of wavelength)
• inﬁnitely thick coating
• diﬃcult to implement constraints on the refractive index
• can not handle the cases when more than one incident angle or polarization is speciﬁed
• can not handle phase speciﬁcations

4.2

Smooth by Penalty Terms

˜
Straightforward optimization of the merit functions R(n, k) − R(k) will lead to solutions
which resemble multilayers or contain many discontinuities. One method to force smoothness
of the solutions is by adding penalty terms. This is the approach adopted in [21], where the
the refractive-index proﬁle n(z) is parameterized as polylines {(zi , ni ) : 1 ≤ i ≤ N } and the
objective function is deﬁned as

1
f (X) =
L

L
j=1

˜
R(X, λj ) − R(λj )
∆j

2

+α

N −1
j=1

nj+1 − nj 2
.
zj+1 − zj

where X = (z1 , · · · , zN , n1 , · · · , nN ). Our approach is similar except n(z) is parameterized
as step functions {(ti , ni ) : 1 ≤ i ≤ N } and use of diﬀerent metrics in the objective functions:

f (X) =


1 L
L

j=1

˜
R(X, kj ) − R(kj )
∆j

1

p p



+α




1
N − 1

N −1
j=1

2(nj+1 − nj )
ti + ti+1

1

q q

.



˜
where p, q are even natural numbers (p can be odd if R(kj ) ≡ 0), ∆j are design tolerances
and α is a weight factor controlling balance between the performance and smoothness of the
36

1.0
2.25
0.8
1.45
0.6
2.25

Refractive Index

0.4
1.45
2.25
0.2

Reﬂected Power

1.0

0.5
0.0

0.8
0.6

0.5
0.0

0.4
0.2

1.45
0.0
0.0
0
2
7
9
0.0 1 0.2 3 4 5 60.6 8 0.8 10 1.04
0.4
0.0
Thickness

0.25

0.4 6 0.6 70.8
Wavelength

0.5
0.0
8
1.0

Figure 4.2: Refractive-index proﬁles (left) and reﬂectance curves (right) of optimal solutions
for a ramp ﬁlter. Top to bottom: α = 0, 0.01, 0.1.
resultant refractive-index proﬁles. Notice that the penalty term is a discretization of the
Lq norm of the n (z). The gradient of f (X) can be computed easily using the derivative
formulas for R derived in chapter 1.
My approach diﬀers from [21] in three aspects. First, n(z) is parameterized by step
functions instead of polylines. Using step functions simpliﬁes the calculation and the results
are similar to those using polylines. Second, eﬀect of diﬀerent metrics in both the merit
(controlled by p) and penalty functions (controlled by q) will be investigated. In particular,
the use of high values for q yield interesting solutions never seen before. Third, the penalty
term is normalized with respect to the number of parameters so that the choice of the value
for α will be independent of the level of discretization.
Figure 4.2 shows the refractive-index proﬁles (left) and reﬂectance curves (right) of optimal solutions for a ramp ﬁlter with diﬀerent values of α. We have set p = q = 2 in this test.
Observe how the refractive-index proﬁle becomes smoother as α increases. As expected, the
performance deteriorates as α increases. Remarkably, a small value of α can smooth the
solution eﬀectively with little eﬀect on the performance.
The eﬀect of p on the solutions are similar to multilayer cases as explained in the previous
37

Reﬂected Power

Refractive index
2.25
1.45

0 1 2 3 4 5 6 7 8 9 10
Thickness

4

5

6
7
Wavelength

8

1.0
0.5
0.0

Figure 4.3: Refractive-index proﬁle (left) and reﬂectance curve (right) of the optimal solution
for a ramp ﬁlter with α = 0.01, p = 2, q = 10.

1.0
2.25
0.8
1.45
0.6
2.25

Refractive Index

1.0

Reﬂected Power
0.5
0.0

0.8
0.6

0.4
0.4
1.45
2.25
0.2
0.2
1.45
0.0
0.0
0
4
18
0.0 2 0.2 6 8 10 12 14 160.8 20 1.04
0.4
0.6
0.0
Thickness

0.5
0.0
0.25

0.4 6 0.6 70.8
Wavelength

0.5
0.0
8
1.0

Figure 4.4: Dependence of the solution on N . Top to bottom: N = 64, 128, 256.
section. The eﬀect of q is more interesting.
Figure 4.3 shows the solutions of the same problem but with α = 0.01 and q = 10. Notice
the solution is qualitatively similar to the q = 2 but the curved parts are straightened up.
On the other hand, the diﬀerence in performance is negligible. Again, this is beneﬁcial since
diﬀerent solutions are always desirable.
Next we investigate the impact of the level of discretization on the performance of the
design. In principle, the ﬁner the discretization the better performance is expected since we
have more degrees of freedom. In practice, however, the impact becomes negligible at rather
rough discretization. Figure 4.4 shows the refractive-index proﬁles and the corresponding
reﬂectance curves for the same problem but with diﬀerent level of discretization.

Notice

the performance of these designs are almost indistinguishable from each other although the
38

refractive-index proﬁles are very diﬀerent. There are two heuristic explanations for this
phenomenon. First, it is a well known fact that small features of a structure has little
eﬀect on waves with large wavelength. Thus reﬁnement of the refractive-index in ﬁne scale
has little impact on its optical performance. This is conﬁrmed from the numerical results,
where we see the refractive-index proﬁles obtained with ﬁner discretization are generally
the interpolations of those obtained with rough discretization. Second, the number of local
minima increases fast as the number of variables increases. So it is more likely to get stuck
in a bad local minimum with ﬁner discretization.
We are motivated by the fact that the refractive-index proﬁles obtained with ﬁner discretization are generally the interpolations of those obtained with rough discretization and
devised a method to speed up the convergence:
1. set N to a small number, e.g. 8
2. optimize until convergence is reached
3. divide each layer into two sub-layers of equal thickness
4. set refractive index of the sub-layers by interpolation from neighbor layers
5. go to step 2 and repeat until desired number of layers is reached
This simple algorithm can increase convergence speed about three fold for most experiments
we conducted. Further improvements are achieved if we increase the convergence tolerance
progressively in step 2.

39

Chapter 5
Global Basis Functions
Existing optimization methods for the optimal design of thin-ﬁlm coatings have used only
local basis functions (step and hat functions) for the parameterization of the refractiveindex proﬁle. Global basis functions such as sinusoidal functions have been used only in
heuristic design methods, where there are only a few design parameters such as the amplitude,
frequency and shifting. The coatings that can be designed by heuristic methods are very
limited.
The Fourier transform method ﬁnds the refractive-index proﬁles as an inverse Fourier
transform of an appropriate function. However, it is not an optimization method and has
many limitations.

5.1

Formulation

The problem of using global basis functions is that there is no direct way to translate the
bound constraints on the refractive-index proﬁle to its spectral coeﬃcients. This is also the
major diﬃculty and limitation of the Fourier transform method.
To overcome this diﬃculty, we transform the set of bounded functions {n(z) ∈ L∞ (0, L) :
nL ≤ n(z) ≤ nH } to the set of unbounded functions {m(z) : m(z) ∈ L∞ (0, L)} through the

40

mapping
m(z) =

w(n(z) − c)

w2 − (n(z) − c)2

,

(5.1)

where w = (nH − nL )/2 and c = (nH + nL )/2. This transformation is obtained by mapping
the line segment (nL , nH ) to (−∞, ∞) through a half circle. This mapping can be inverted
to get an expression of n(z) in terms of m(z):

wm(z)

n(z) = c +

m(z)2 + w2

.

(5.2)

This particular choice of mapping is picked by intuition and is not necessarily the best.
In fact, there are inﬁnitely many choices out there, e.g. the tangent function

m(z) = tan

π
(x − c) .
w

Now there is no bound constraints on m(z), let us approximate it by truncated Fourier
series
N

m(z) ≈

aj cos
j=0

jπ
x +
L

N

bj sin
j=1

jπ
x .
L

(5.3)

It is certainly possible to use other global basis functions such as Chebyshev or Hermite
polynomials.
Given a vector of coeﬃcients {a0 , · · · , aN , b1 , · · · , bN } we have a function m(z) from (5.3)
and a refractive-index proﬁle n(z) by solving (5.1). The reﬂectance curve corresponding to
n(z) can be then computed from the matrix method. Thus we have a function mapping
the Fourier coeﬃcients to reﬂectance and the unconstrained optimization problem for a

41

broadband coating:

min

a0 ,··· ,aN ,b1 ,··· ,bN

˜
R(k; a0 , · · · , aN , b1 , · · · , bN ) − R(k) Lp (k ,k ) ,
L H

(5.4)

where R(k; a0 , · · · , aN , b1 , · · · , bN ) is the reﬂectance at wavenumber k of a refractive-index
˜
proﬁle n(z) obtained from the coeﬃcients {a0 , · · · , aN , b1 , · · · , bN } and R(k) is the target
reﬂectance function.
To compute the reﬂectance and its gradient, we discretize n(z) and m(z) by piece-wise
constant functions. Divide the coating into D sub-layers of equal thickness. Denote zi the
center of the i − th layer and ni = n(zi ), mi = m(zi ), Then the reﬂectance R is computed
from the matrix method and the partial derivatives of R is given by

∂R
=
∂aj

D
i=1
D

=
∂R
=
∂bj

i=1
D
i=1
D

=
i=1

∂R ∂ni ∂mi
∂ni ∂mi ∂aj

(5.5a)

3
∂R 3 2
2 )− 2 cos
w (mi + w
∂ni

jπ
z , 0 ≤ j ≤ N,
L i

∂R ∂ni ∂mi
∂ni ∂mi ∂aj

(5.5b)

(5.5c)

3
∂R 3 2
2 )− 2 sin
w (mi + w
∂ni

jπ
z , 1 ≤ j ≤ N,
L i

(5.5d)

where ∂R/∂ni is also computed from the matrix method.

5.2

Numerical Experiments

Consider a broadband antireﬂection coating over the wavelength range [λL , λH ] = [400, 1200]
at TE polarization and normal incident angle. The refractive indices of the incident medium
42

1.0
2.35
0.8
1.45
2.35
0.6
0.4
1.45
2.35
0.2
1.45
0.0
0.0
0

Refractive Index

1.0

Reﬂected Power

0.8
0.6
0.4
0.2
0.2
200

0.0
0.4
0.0
400 0.6
600 0.8 800 1.0
400
Thickness

0.2
0.8
600 0.4 800 0.6 1000
Wavelength

0.02
0.01
0.00
0.02
0.01
0.00
0.02
0.01
0.00
1.0

Figure 5.1: Snapshots of the refractive-index proﬁle (left) and reﬂectance (right) at 10, 50,
and 200 iterations (top to bottom) during the optimization of a broadband AR coating
parameterized by sinusoidal functions.
and substrate are na = 1.0 and ns = 1.52 respectively. The lower and upper bound of the
refractive index of the coating is nL = 1.45 and nH = 2.35 respectively. The total thickness
of the coating is L = 900. The number of terms in the Fourier series is N = 20 and the
initial guess is {a0 , · · · , aN , b1 , · · · , bN } = {1, · · · , 0, 0, · · · , 0}, corresponding to m(z) = 1
and n(z) = (nL + nH )/2. The Matlab function “fmincon” is used to run the optimization.
Figure 5.1 shows the snapshots of the refractive-index proﬁle n(z) and the corresponding
reﬂectance curves R(λ) at 10, 50, and 200 iteration steps during the optimization.
We have a number of observations from this example.
First, it is interesting to see how n(z) evolves from a smooth function to an almoststep function taking only two values, namely the lower and upper bounds nL , nH . This is
consistent with the maximum principle of thin-ﬁlm optics [38].
Second, n(z) varies much more signiﬁcant than R(λ) does during the iterations. This is
due to the ill-posedness of the corresponding inverse problem of recovering n(z) from R(λ).
While ill-posedness is generally considered to be a disadvantage of inverse problems, it is
the opposite in this case since we can have very diﬀerent designs with similar performance.

43

For example, the rugate (2nd) and multilayer (3rd) designs in Figure 5.1 are qualitatively
diﬀerent but their performance are almost identical.
Third, the reﬂectance curves are uniformly distributed in the spectral region even if p = 1
is used in the objective function (5.4). In contrast, the reﬂectance curves obtained by using
local basis functions with small values of p possess a peak at the upper bound of the spectral
region.
Ongoing and future works in this topic include the possibilities to use other transformations and basis functions, determine the appropriate value for the number of terms N or
how to increase it progressively to obtain better solutions or increase convergence speed. We
shall also consider the use of penalty or constraint functions to force smoothness of n(z) as
we did in chapter 4.

44

Chapter 6
Diﬀraction Gratings
So far we considered only optical coatings built from layer structures. The refractive index
of such structures vary along one direction and is homogeneous over the other two. By
allowing the refractive index of the structure to vary periodically along one or both of
the other two directions one obtains singly or doubly periodic structures, often called 2D
or 3D diﬀraction gratings, respectively. Dobson studied the optimal design of wide-angle
antireﬂective diﬀraction gratings [7] and blazed diﬀraction gratings [8]. Following the spirit
of these work, we study the optimal design of broadband optical coatings by the use of
diﬀraction gratings. The value and gradient of the objective function will be obtained by
solving a 2D PDE with nonlocal transparent boundary conditions [9] and the corresponding
adjoint equation [7], respectively.

6.1

Model and Forward Solver

Consider a 1D diﬀraction grating so that µ = 1 and ε is invariant along y-axis with



ε ,
 1
z ≥ b − δ,





ε(x, z) = ε(x + Λ, z), −b + δ < z < b − δ,







ε 2 ,

z ≤ −b + δ

45

z
y
x

θ

n(x, z) = n(x + Λ, z)

···

···

Λ
Figure 6.1: Problem geometry for 1D diﬀraction gratings.
where b > 0, δ > 0 are a constants with b − δ > 0. Figure 6.1 summarizes the geometry for
1D diﬀraction gratings.
Given a time harmonic plane incident wave, we will determine the reﬂectance and transmittance by the variational approach with transparent boundary conditions [9]. The Maxwell’s
equations for a time harmonic wave with frequency ω reduce to the following coupled system
of linear equations

× E − ikH = 0,

(6.1a)

× H + ikεE = 0,

(6.1b)

where k = ω/c is the wavenumber in vacuum. We also have the following boundary conditions at a surface of discontinuity:

n12 × (E2 − E1 ) = 0,

(6.2a)

n12 × (H2 − H1 ) = 0,

(6.2b)

46

In the case of TE polarization, the system of linear equations reduces to the Helmholtz
equation
(∆ + k 2 ε)u = 0,
where u is the tangential component of E. For the sake of uniqueness, we only look for quasi
periodic solutions u such that
uα = u exp(−iαx)
is Λ periodic, where α = kε1 sin θ and θ is the angle of incidence.
Denote

Ω = {(x, z) : −b < z < b},

(6.3a)

Ω1 = {(x, z) : z ≥ b},

(6.3b)

Ω2 = {(x, z) : z ≤ −b},

(6.3c)

Γ1 = {(x, z) : z = b},

(6.3d)

Γ2 = {(x, z) : z = −b}.

(6.3e)
(6.3f)

For each n ∈ Z denote αn = 2πn/Λ. By expanding uα into its Fourier series along the
x-axis:

(n)

uα (x, z) =

uα (z) exp(iαn x),

(6.4a)

n∈Z
(n)

uα (z) =

1 Λ
uα (x, z) exp(−iαn x)dx,
Λ 0

47

(6.4b)

the fact that uα is a sum of plane waves in Ω1 , Ω2 :

uα (x, z)|Ωj =

(n)

aj

(n)

exp[±iβj z + iαn x],

(6.5)

n∈Z

and the continuity of u and ∂u/∂ν across Γ1 , Γ2 , we obtain the following boundary value
problem for uα

(∆α + k 2 ε)u = 0, in Ω,

(6.6a)

∂
− T1 )uα = −2iβ exp(−iβb), on Γ1 ,
∂ν
∂
− T2 )uα = 0, on Γ2 ,
(
∂ν
(

(6.6b)
(6.6c)

where
∆α = ∆ + 2iα∂x − |α|2 ,
(n)

iβj f (n) exp(iαn x),

Tj f =
n∈Z

f (n) =
(n)

βj

=

1 Λ
f (x) exp(−iαn x)dx,
Λ 0
k 2 εj − (αn + α)2

(6.7)
(n)

s.t. arg(βj ) ∈ [0, π),

β = kε1 cos θ
Here and in the sequel the domain and boundary will be restricted into one periodic cell
{(x, z) : 0 < x < Λ}.
Denote H(Ω) = {v : v ∈ H 1 (Ω), v(0) = v(L)} where H 1 (Ω) is the standard Sobolev
space. Multiplying the PDE in (6.6a) by the conjugate of a test function φ ∈ H(Ω) and
applying Green’s formula with the given boundary conditions, we obtain the variational, or

48

weak, formulation of the boundary value problem (6.6a):

−

Ω

+
Γ1

uα ·

φ+

(T1 uα )φ +

Ω
Γ2

(k 2 ε − α2 )uα φ + 2iα

Ω

(∂x uα )φ

(T2 uα )φ

exp(−iβb)φ.

= 2iβ
Γ1

This variational formulation will be discretized and solved by the ﬁnite element method.
The layered structures studied in the previous chapters have the property that the reﬂected light propagates in a single direction, so is the transmitted light. On the other hand,
diﬀraction gratings can have light reﬂected and transmitted in multiple directions. This can
(n)

be seen from (6.5) and the deﬁnition of βj

(n)

in (6.7). Each real value of βj

corresponds to

a propagating mode and each complex value corresponds to an evanescent mode (propagates
along the x-axis and decreases exponentially along z-axis).
(n)

For real values of j , j = 1, 2, the set of values of n such that βj

Λj = {n ∈ Z : |αn + α| < k 2 εj }.

is real is given by

(6.8)

Once the solution uα is obtained, we can determine for each order n ∈ Λj the amplitude
reﬂectance rn and amplitude transmittance tn from (6.4a) and (6.5):

(0)

(0)

(n)

(n)

(0)

r0 = uα (b) exp(−iβ1 b) − exp(−2iβ1 b),
rn = uα (b) exp(−iβ1 b),
(n)

(n)

tn = uα (−b) exp(−iβ2 b),

49

n ∈ Λ1 ,
n ∈ Λ2 .

(6.9a)
(6.9b)
(6.9c)

As for the layered structures, we then use Poynting vector to determine the corresponding
power reﬂectance Rn and power transmittance Tn and obtain the following:

Rn = |rn
Tn = |tn

(n)
2 cos θ1 ,
|

(6.10a)

cos θ

(n)
2 cos θ2
|

ε2
,
ε1

cos θ

(6.10b)
(6.10c)

(n)

where θj , j = 1, 2 is the angle of propagation for mode of order n and is given by

(n)

θj

6.2

= arctan 



αn + α 
(n)

.

βj

Adjoint Problem and Gradient

As in the layered structures, we need to ﬁnd a method to compute the gradient of the
reﬂectance and transmittance if we would like to use gradient based optimization algorithms.
Unlike the layered structures, we do not have a explicit expression relating the reﬂectance
Rn to the permittivity ε. We derive the formula for the gradient by extending the formal
argument as in [7].
Let ε be ﬁxed and δε be a small perturbation. By (6.4a) and (6.9a) we obtain the
directional derivative

Dε rn (δε) =

(n)
2 exp(−iβ1 b)
k

Λ

Γ1

50

δuα exp(−iαn x)dx,

(6.11)

where δuα is the solution for the following boundary value problem:

(∆α + k 2 ε)δuα = −k 2 δεuα ,
Tj −

∂
∂ν

δuα |Γj = 0,

j = 1, 2.

(6.12a)
(6.12b)

Let wα be the solution of the following “adjoint problem”:

(∆α + k 2 ε)wα = 0,

(6.13a)

∗
T1 −

∂
∂ν

wα |Γ1 = −φ(x),

(6.13b)

∗
T2 −

∂
∂ν

wα |Γ2 = 0,

(6.13c)
(6.13d)

where φ(x) = rn exp(iαn x) and

(n)

∗
Tj f =
n∈N

−iβj f (n) exp(iαn x).

Multiplying the ﬁrst equation in (6.12a) by wα and applying Green’s formula yields

Ω

δεuα wα =

Γ1

= rn
=

δuα φ

Γ1

δuα exp(−iαn x)
Λ
(n)

k 2 exp(−iβ1 )b

(6.14)

rn Dε rn (δε).

Hence
Dε Rn (δε) =

(n)
2 exp(−iβ1 )b
2k

Λ

51

Ω

δεuα wα .

(6.15)

This implies the L2 -gradient of Rn with respect to ε is given by
(n)

(n)

ε Rn

=

The expressions for

6.3

2k 2

exp(iβ1 )b
exp(−iβ1 )b
uα wα = 2k 2
uα wα ,
Λ
Λ
Ω
Ω

ε Tn

(6.16)

can be obtained in a similar fashion.

Numerical Experiments

In order to solve the optimization problem numerically, we need ﬁrst discretize the functional
space L∞ (Ω) for ε. We are going to use ﬁnite element method with rectangular bilinear
elements, so it is most natural here to discretize ε by piece-wise constant functions over
rectangular mesh. For real ε we deﬁne

εi,j = χ[(j−1)∆x ,j∆x ]×[(i−1)∆y ,i∆y ] ,

where χ is the standard characteristic function. For complex valued ε we can consider the
real and imaginary parts separately. For the sake of convenience, the mesh size ∆x , ∆y are
chosen to be integer multiples of the mesh size for the ﬁnite element space. By (6.16) we
obtained the discrete version of the gradient

i,j Rn

=



(n)


exp(iβ1 )b
2k 2
uα wα εi,j .


Λ
Ω

(6.17)

Let us consider the design of a diﬀraction grating with broadband antireﬂection proper-

52

ties. The objective function is a normalized discretization of R0 1 :
1
f( ) =
N

where

N

R( ; ki ),
i=1

= {εi,j } are design variables and {ki } is a uniform set of sampling points in the

spectral region of interest.
The refractive indices of the incident medium and substrate are 1.0 and 1.52 respectively. The lower and upper bound of the refractive index of the coating are 1.45 and 2.35
respectively. We assume normal incidence and TE polarization. The range of wavelength is
√
[λL , λH ] = [1, 3]. The depth of the grating is 2b = 2 and the period Λ = 0.99λL ε1 . Notice
this value of Λ will guarantee only the zero-th reﬂecting order is propagating so that only R0
is nonzero for all wavelength by (6.8). By increasing Λ we will introduce more propagating
orders and it is up to the user to minimize one of the orders or the combination of them.
The domain Ω is divided into 32 × 64 uniform rectangular mesh and ε is discretized by
16 × 32 uniform rectangular mesh. The minimization of the objective function is done by
the “fmincon” command in Matlab, and the algorithm is chosen to be “active-set”. The
maximal iteration step is set to 300.
In the ﬁrst experiment we choose the initial guess to be the homogeneous proﬁle ε =
√

ε1 ε2 . Figure 6.2 shows the refractive-index proﬁles and reﬂectance curves for the initial

guess and optimized solution. The refractive index is normalized from [nL , nH ] to [0, 64]
so that it can be represented as a gray scale image. This transformation is used in all the
following ﬁgures.
A notable feature of the optimized proﬁle it resembles a multilayer structure, i.e the
refractive index in invariant along the x-axis and is mostly piece-wise constants along the

53

1
60
5

0.8

50
40
15
30
20

Reflectance

z−axis

10
0.6
0.4

20
25

0.2

10
30
5

10
x−axis

0
1

15

1.5

2
2.5
Wavelength

3

1.5

2
2.5
Wavelength

3

0.2
60
5
50

0.15

40
15
30
20

Reflectance

z−axis

10

20

0.1

0.05

25
10
30
5

10
x−axis

0
1

15

Figure 6.2: Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection
coating with the homogeneous proﬁle as the initial guess. Top: initial guess; Bottom: optimized.

54

z-axis. This is expected and proofs the validity of our computation of the gradient and the
algorithm. The initial guess is a layered structure, therefore the gradient of R0 should be
the same as if we are solving a layered problem. So the gradient will be zero along the x-axis
and the next iteration will stay as a layered structure. This implies we can never achieve
better performance if the initial guess is a multilayer structure. Indeed, the reﬂectance curve
in Figure 6.2 is similar to and no better than those obtained from the optimization of a pure
multilayer coating with similar physical parameters.
Therefore to obtain a solution that is not a multilayer, we need select the initial guess to
be non-multilayer. In the next example, we choose the initial guess to be have a rectangular
proﬁle, i.e.
ε(x, z) =




ε ,
 L

x ∈ (0, Λ/2),



ε , x ∈ (Λ/2, Λ).
 H
The refractive-index proﬁles and reﬂectance curves for the initial guess and optimized
solution are shown in Figure 6.3. First of all, notice the distinct feature of the reﬂectance
curve for the initial guess, i.e. the sharp peaks. The existence of these peaks is due to
the resonance eﬀect of diﬀraction gratings. The interested reader is referred to Bao and
Ramdani [18], Bao and Huang [16, 17] and the references cited therein for the theoretical
justiﬁcation and applications of this eﬀect. The references just mentioned also considered
the optimal design of the resonance eﬀect, i.e. ﬁnd a grating proﬁle so that the peaks appear
at desired locations and whose width are as small as possible.
Notice the peaks are not present in the reﬂectance curves when the proﬁle is homogeneous.
This is because a diﬀraction grating with a homogeneous proﬁle is indeed not a periodic but
layered structure, so the resonance eﬀect will not occur.

55

1
60
5

0.8

50
40
15
30
20

Reflectance

z−axis

10
0.6
0.4

20
25

0.2

10
30
5

10
x−axis

0
1

15

1.5

2
2.5
Wavelength

3

1.5

2
2.5
Wavelength

3

1
60
5

0.8

50
40
15
30
20

Reflectance

z−axis

10
0.6
0.4

20
25

0.2

10
30
5

10
x−axis

0
1

15

Figure 6.3: Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection
coating with the rectangular proﬁle as the initial guess. Top: initial guess; Bottom: optimized.

56

For the purpose of antireﬂection coatings, the peaks should be suppressed as much as
possible. This is indeed the case in the for the optimized proﬁle in Figure 6.3, where not only
the peaks are suppressed but also the average reﬂectance is reduced signiﬁcantly. In fact,
the reﬂectance curve is pretty much on the same level as that obtained with a homogeneous
initial guess. On the other hand, the refractive-index proﬁle is not a multilayer and very
diﬀerent from that obtained with a homogeneous initial guess.
Next we choose the initial guess to be a triangular proﬁle, i.e.

ε(x, z) =




ε ,
 L


ε ,
 H

z>

4b
Λ

x−

Λ
4

or z > −

4b
Λ

x−

3Λ
4

,

else.

Figure 6.4 shows the refractive-index proﬁles and reﬂectance curves for the initial guess and
optimized solution. Like the case when the initial guess is rectangular, the reﬂectance curve
in this case also possesses sharp peaks. They are suppressed in the optimized solution, but
not as well as in the previous experiment. The average reﬂectance, however, is pretty much
the same as other examples.

57

1
60
5

0.8

50
40
15
30
20

Reflectance

z−axis

10
0.6
0.4

20
25

0.2

10
30
5

10
x−axis

0
1

15

1.5

2
2.5
Wavelength

3

1.5

2
2.5
Wavelength

3

1
60
5

0.8

50
40
15
30
20

Reflectance

z−axis

10
0.6
0.4

20
25

0.2

10
30
5

10
x−axis

0
1

15

Figure 6.4: Refractive-index proﬁles and reﬂectance curves for a broadband antireﬂection
coating with the triangular proﬁle as the initial guess. Top: initial guess; Bottom: optimized.

58

Chapter 7
Future Work
Much of the work in this dissertation is not completely ﬁnished and there are quite a few
open problems and new directions that I would like to consider in the future. The topics
include but are not limited to the following
• Theoretical justiﬁcation of the observations made on the eﬀect of the objective functions. We shall investigate in depth the properties of the objective functions in the
context of optical coatings. We hope the investigations not only help us understand the
eﬀect but also provide guidelines and even new algorithms for numerical computations.
• Extend the study on the eﬀect of objective functions to other types of coatings such as
mirrors, bandpass ﬁlters and edge ﬁlters. We believe similar eﬀects will be observed.
• Use of diﬀerent metrics in diﬀerent parts of the spectral or angular region of interest.
In the experiments with broadband AR coatings, for example, the reﬂectance is usually
larger towards the longer wavelength region if p is small. By using larger value for p,
the reﬂectance is more evenly spread out, at the cost of higher reﬂectance in the shorter
wavelength region and a higher average value. We propose using smaller values for p
in the shorter wavelength region and larger values in the longer wavelength region. We
hope this will suppress the heightening and keep the average low at the same time.
• In the design of rugate coatings by penalty terms, I would like to derive theoretical
59

or empirical expressions relating the performance to the parameters, especially the
weight factor α. Such expressions enable users to determine the best value of α for
their practical needs.
• In the design of rugate coatings by penalty terms, we found the optimization algorithm
converges faster as the weight factor α increases. Since the convergence is faster for
larger α, I also plan to devise algorithms in which α is gradually decreased during the
optimization in order to increase the speed of convergence for small α or ﬁnd better
solutions.
• We have already seen that good rugate coatings can be obtained by stopping the
iteration prematurely when we use global basis functions to discretize the refractiveindex proﬁle. We shall consider adding a penalty term to see what solutions we can
obtain as in the case where the proﬁle is discretized by local basis functions.
• The diﬀraction gratings considered in this dissertation have the most general proﬁle,
i.e a 2D function. This poses diﬃculties in terms of both feasibility and numerical computations. It seems more practical to consider diﬀraction gratings deﬁned by interface,
i.e. a 1D function. The optimal design of such type of gratings have been considered
in [15, 8, 14, 13] but only at a single frequency. We shall consider the broadband,
omnidirectional and other more complicated situations.
• Optimal design problems are closely related to inverse problems. Recently a study [36]
showed the methods for inverse problems can be successfully applied for the optimal
design of layered structures. I will consider extending them to the optimal design of
diﬀraction gratings.

60

• Recently it has been shown the methods of recursive linearization are very eﬀective
for inverse problems in diﬀractive optics [39, 40]. The problems considered there use
input data at multiple frequencies, thus it is similar to the optimal design of broadband
coatings. I plan to adopt this idea to the optimal design of not only diﬀraction gratings
but also layered structures.
• Optimal design of 3D diﬀraction gratings. The computational cost will increase dramatically and we shall investigate numerical tricks such as increasing the discretization
level and convergence tolerance progressively.

61

BIBLIOGRAPHY

62

BIBLIOGRAPHY
[1] Sh. A. Furman and A. V. Tikhonravov. Basics of Optics of Multilayer Systems. GifSur-Yvette, Edition Frontieres, 1992.
[2] H. A. Macleod. Thin-Film Optical Filters, Third Edition. Taylor & Francis, 2001.
[3] P. Baumeister. Design of multilayer ﬁlters by successive approximations. J. Opt. Soc.
Am., 48(12):955–957, Dec 1958.
[4] J. A. Dobrowolski and R. A. Kemp. Reﬁnement of optical multilayer systems with
diﬀerent optimization procedures. Appl. Opt., 29(19):2876–2893, Jul 1990.
[5] A. V. Tikhonravov. On the synthesis of optical thin ﬁlms using optimality conditions.
Vestn. Mosk. Univ. Fiz. Astronomiya, 23:91–93, 1982.
[6] A. V. Tikhonravov. Some theoretical aspects of thin-ﬁlm optics and their applications.
Appl. Opt., 32(28):5417–5426, Oct 1993.
[7] D. C. Dobson. Optimal design of periodic antireﬂective structures for the helmholtz
equation. European Journal of Applied Mathematics, 4:321–339, 1993.
[8] D. C. Dobson. Optimal shape design of blazed diﬀraction gratings. Appl. Math. Optim.,
40:61–78, 1999.
[9] G. Bao, D. C. Dobson, and J. A. Cox. Mathematical studies in rigorous grating theory.
J. Opt. Soc. Am. A, 12(5):1029–1042, May 1995.
[10] M. Born and E. Wolf. Principles of Optics: Electromagnetic Theory of Propagation,
Interference and Diﬀraction of Light (7th Edition). Cambridge University Press, 7th
edition, 1999.
[11] P. G. Verly. Modiﬁed needle method with simultaneous thickness and refractive-index
reﬁnement for the synthesis of inhomogeneous and multilayer optical thin ﬁlms. Appl.
Opt., 40(31):5718–5725, Nov 2001.
[12] C. Sanderson. Armadillo: C++ linear algebra library. http://arma.sourceforge.net.

63

[13] G. Bao. Inverse and optimal design problems in diﬀractive optics. In Y. Hon, M. Yamamoto, J. Cheng, and J. Lee, editors, Recent Developments in Theories & Numerics,
Proceedings of the International Conf. on Inverse Problems, pages 7–46. World Sci.,
Hong Kong, 2003.
[14] G. Bao and E. Boneetier. Optimal design of periodic diﬀractive structures. Appl. Math.
Optim., 43:103–116, 2001.
[15] G. Bao and D. C. Dobson. Modeling and optimal design of diﬀractive optical structures.
Survey on Industrial Math, 8:37–62, 1998.
[16] G. Bao and K. Huang. Optimal design of guided mode grating resonance ﬁlters. IEEE
Photonic Tech. Lett., 16(1):141–143, 2004.
[17] G. Bao and K. Huang. Computational design for guided-mode grating resonances. J.
Opt. Soc. Am. A, 22(7):1408–1413, Jul 2005.
[18] G. Bao and K. Ramdani. Resonant frequencies for diﬀraction gratings. Applied Mathematics Letters, 15:755–760, 2002.
[19] G. Bao, K. Huang, and G. Schmidt. Optimal design of nonlinear diﬀraction gratings.
J. Comput. Phys., 184:106–121, 2003.
[20] J. A. Dobrowolski. Completely automatic synthesis of optical thin ﬁlm systems. Appl.
Opt., 4(8):937–946, Aug 1965.
[21] A. V. Tikhonravov, M. K. Trubetskov, T. V. Amotchkina, M. A. Kokarev, N. Kaiser,
O. Stenzel, S. Wilbrandt, and D. G¨bler. New optimization algorithm for the synthesis
a
of rugate optical coatings. Appl. Opt., 45(7):1515–1524, March 2006.
[22] R. Fletcher. Practical methods of optimization. Wiley, New York, 2nd edition, 1987.
[23] G. Erdmann and F. Santosa. Minimax design of optically transparent and reﬂective
coatings. J. Opt. Soc. Am. A, 21(9):1730–1739, Sep 2004.
[24] A. Premoli and M. L. Rastello. Minimax reﬁning of optical multilayer systems. Appl.
Opt., 31(10):1597–1605, Apr 1992.
[25] A. Premoli and M. L. Rastello. Minimax reﬁning of wideband antireﬂection coatings
for wide angular incidence. Appl. Opt., 33(10):2018–2024, Apr 1994.

64

[26] S. G. Johnson. The nlopt nonlinear-optimization package. http://ab-initio.mit.
edu/nlopt.
[27] D. C. Liu and J. Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical Programming, 45:503–528, 1989. 10.1007/BF01589116.
[28] D. Kraft. Algorithm 733: Tompfortran modules for optimal control calculations. ACM
Transactions on Mathematical Software, 20(3):262–281, 1994.
[29] P. G. Verly. Optical coating synthesis by simultaneous refractive-index and thickness
reﬁnement of inhomogeneous ﬁlms. Appl. Opt., 37(31):7327–7333, Nov 1998.
[30] T. V. Amotchkina. Empirical expression for the minimum residual reﬂectance of normaland oblique-incidence antireﬂection coatings. Appl. Opt., 47(17):3109–3113, 2008.
[31] A. V. Tikhonravov, M. K. Trubetskov, T. V. Amotchkina, and J. A. Dobrowolski.
Estimation of the average residual reﬂectance of broadband antireﬂection coatings. Appl.
Opt., 47(13):C124–C130, 2008.
[32] R. Willey. Predicting achievable design performance of broad-band antireﬂection coatings. Appl. Opt., 32:5447–5451, 1993.
[33] R. Willey. Reﬁned criteria for estimating limits of broad-band ar coatings. Proc. SPIE,
5250:393–399, 2004.
[34] J. A. Dobrowolski, A. V. Tikhonravov, M. K. Trubetskov, B. T. Sullivan, and P. G. Verly.
Optimal single-band normal-incidence antireﬂection coatings. Appl. Opt., 35(4):644–
658, Feb 1996.
[35] A. V. Tikhonravov and M. K. Trubetskov. Modern design tools and a new paradigm in
optical coating design. Appl. Opt., 51(30):7319–7332, Oct 2012.
[36] S. W. Anzengruber, E. Klann, R. Ramlau, and D. Tonova. Numerical methods for the
design of gradient-index optical coatings. Appl. Opt., 51(34):8277–8295, Dec 2012.
[37] B. G. Bovard. Rugate ﬁlter theory: an overview. Appl. Opt., 32(28):5427–5442, October
1993.
[38] A. V. Tikhonravov, M. K. Trubetskov, and G. W. DeBell. Application of the needle
optimization technique to the design of optical coatings. Appl. Opt., 35(28):5493–5508,
Oct 1996.
65

[39] G. Bao, P. Li, and J. Lv. Numerical solution of an inverse diﬀraction grating problem
from phaseless data. J. Opt. Soc. Am. A, 30(3):293–299, March 2013.
[40] G. Bao, P. Li, and H. Wu. A computational inverse diﬀraction grating problem. J. Opt.
Soc. Am. A, 29(4):394–399, April 2012.

66