HAL Id: hal-00867708
https://hal-iogs.archives-ouvertes.fr/hal-00867708
Submitted on 30 Sep 2013
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Numerical performance of nite-dierence modal
methods for the electromagnetic analysis of
one-dimensional lamellar gratings
Philippe Lalanne, Jean-Paul Hugonin
To cite this version:
Philippe Lalanne, Jean-Paul Hugonin. Numerical performance of nite-dierence modal methods for
the electromagnetic analysis of one-dimensional lamellar gratings. Journal of the Optical Society of
America. A Optics, Image Science, and Vision, Optical Society of America, 2000, 17 (6), pp.1033-1042.
�hal-00867708�
Numerical performance of finite-difference modal
methods for the electromagnetic
analysis of one-dimensional lamellar gratings
Philippe Lalanne and Jean-Paul Hugonin
Laboratoire Charles Fabry de l’Institut d’Optique The
´
orique et Applique
´
e, Centre National de la Recherche
Scientifique, B.P. 147, 91403 Orsay Cedex, France
Received October 7, 1999; accepted January 10, 2000; revised manuscript received February 16, 2000
The numerical performance of a finite-difference modal method for the analysis of one-dimensional lamellar
gratings in a classical mounting is studied. The method is simple and relies on first-order finite difference in
the grating to solve the Maxwell differential equations. The finite-difference scheme incorporates three fea-
tures that accelerate the convergence performance of the method: (1) The discrete permittivity is interpolated
at the lamellar boundaries, (2) mesh points are located on the permittivity discontinuities, and (3) a nonuni-
form sampling with increased resolution is performed near the discontinuities. Although the performance
achieved with the present method remains inferior to that achieved with up-to-date grating theories such as
rigorous coupled-wave analysis with adaptive spatial resolution, it is found that the present method offers
rather good performance for metallic gratings operating in the visible and near-infrared regions of the spec-
trum, especially for TM polarization. © 2000 Optical Society of America [S0740-3232(00)00606-2]
OCIS codes: 050.2770, 050.1950, 050.1970.
1. INTRODUCTION
Many rigorous methods exist for analyzing the diffraction
by surface-relief gratings, and much work has been done
during recent years to improve and generalize rigorous
methods. The differential method
1
was associated with
the R-matrix algorithm to improve its stability
2
for highly
conductive metallic gratings, especially for TM polariza-
tion (magnetic field vector perpendicular to the grating
vector); the integral method
3
was generalized to study dif-
fraction by echelles covered with dielectric layers.
4
The
rigorous coupled-wave analysis
5
(RCWA) was improved in
a similar manner,
6
and recently its convergence for TM
polarization and conical mounts was substantially
improved.
7,8
The method of coordinate transformation
9
was extended
10
to gratings with vertical facets. Finite-
difference (FD) methods
11
for solving partial differential
equations are also widely used in electromagnetism for
solving Maxwell’s equations. The finite-difference time
domain (FDTD) method
12
is one example used exten-
sively. FD or finite-element methods are not widely used
in grating theory
13–16
but are often used to study the dif-
fraction by aperiodic objects of finite dimension
17–19
be-
cause of their suitability for incorporating absorbing
boundary conditions to limit the computational domain.
This paper is devoted to a very simple method for the
analysis of one-dimensional (1D) lamellar gratings under
classical mounts. The method shares many features
with standard grating theories such as the differential
method and the RCWA but uses a FD approach to solve
Maxwell’s equations in the grating region. It is twofold.
First, by the application of a FD technique, the modes in-
side the grating region are computed as eigenvectors of a
propagation operator. Then to compute the diffraction
efficiencies, the boundary conditions at the grating inter-
faces are matched by use of the method of moments.
Note that this approach differs strongly from standard
FD methods
13–16,18
that do not require solution of any
eigenproblem but do require a two-dimensional (2D) mesh
for the discretization of a 1D grating. The present
method, in contrast, is time harmonic and uses a 1D mesh
for the discretization of 1D gratings. It is therefore simi-
lar to numerical techniques that are based on FD modal
approaches and used in waveguide theories.
20–23
The FD modal method is described in Section 2 with
special attention devoted to the TM polarization case.
The TE polarization case is briefly reported. In Section 3
the convergence rate of the present method is studied for
metallic and dielectric lamellar gratings for TE and TM
polarizations. It is compared with that achieved by the
classical RCWA and by the enhanced version of the
RCWA recently reported.
24
Limitations of the present
method for analyzing dielectric gratings are emphasized.
Section 4 contains several comments on the present
method. Section 5 concludes the paper.
2. FINITE-DIFFERENCE MODAL METHOD
Let us consider a 1D lamellar grating along the x axis
with a relative permittivity profile
⑀
(x); see Fig. 1. The z
axis is perpendicular to the grating boundaries. The dif-
fraction problem is invariant in the y direction. Mag-
netic effects are not considered in this paper, and the con-
stant
0
denotes the permeability of the periodic
structure.
⑀
0
is the permittivity of the vacuum. The
grating period is denoted by ⌳, and the modulus K of the
grating vector is equal to 2
/⌳. An incident plane wave
of frequency
and wavelength in the vacuum makes an
angle
with the z direction in a nonconical mounting.
P. Lalanne and J.-P. Hugonin Vol. 17, No. 6/June 2000/J. Opt. Soc. Am. A 1033
0740-3232/2000/061033-10$15.00 © 2000 Optical Society of America
We denote the modulus of the wave vector of the incident
wave by k
0
(k
0
⫽ 2
/) in the vacuum. A temporal de-
pendence in e
j
t
of the wave is assumed. In the follow-
ing, the FD modal method is first described for TM polar-
ization. The TE polarization case is then briefly
outlined.
A. TM Polarization
The incident normalized magnetic field is given by
⌿
inc
⫽ exp
兵
⫺jk
0
n
1
关
sin
共
兲
x ⫹ cos
共
兲
z
兴
其
. (1)
The Rayleigh expansions for the magnetic field in the in-
cident medium and in the substrate are given by
⌿
1
⫽ ⌿
inc
⫹
兺
i
R
i
exp
关
⫺j
共
k
xi
x ⫺ k
1,zi
z
兲
兴
, (2)
⌿
3
⫽
兺
i
T
i
exp
兵
⫺j
关
k
xi
x ⫹ k
3,zi
共
z ⫺ h
兲
兴
其
, (3)
respectively. In Eqs. (2) and (3), R
i
and T
i
are the
backward- and forward-diffracted amplitudes, k
xi
is equal
to k
0
n
1
sin
⫺ iK, and k
p,zi
( p ⫽ 1 or 3) is defined by
k
xi
2
⫹ k
p,zi
2
⫽ n
p
2
k
0
2
, with Re(k
p,zi
) ⫹ Im(k
p,zi
)
⬎
0. The
magnetic field ⌿(x) in the grating region satisfies the
Helmholtz equation
2
⌿
z
2
⫹
⑀
x
冉
1
⑀
⌿
x
冊
⫹ k
0
2
⑀
⌿ ⫽ 0. (4)
To solve this differential equation, we use a FD scheme in
the x direction and then calculate analytically the propa-
gation in the z direction. The first step in applying a FD
method is to select a discrete set of values of x (the dis-
crete points) inside one grating period, the x
i
’s, i
⫽ 1,... N, represented with crosses in Fig. 2. The mag-
netic field becomes a function of the discrete index i and
can be noted as a vector ⌿ with N components, the ⌿
i
’s in
Fig. 2. The FD expression for the first derivative along
the x direction is computed at the discrete locations x
i
⬘
.
The x
i
⬘
’s are represented by circles in Fig. 2. An inter-
laced grid similar to the Yee’s space lattice
25
used in
FDTD methods is thus considered. In its discretized
form, Eq. (4) is written as
2
⌿
z
2
⫹ E⌿ ⫽ 0, (5)
where
E ⫽
⑀
1
共
D
2
⑀
2
⫺1
D
1
⫹ k
0
2
I
兲
. (6)
In Eq. (6) the N ⫻ N matrix E is tridiagonal with nonnull
upper right and lower left coefficients, I is the identity
matrix, and
⑀
1
and
⑀
2
are two diagonal matrices that re-
sult from the relative permittivity profile. Because this
profile is piecewise constant, a naı
¨
ve discretization of the
relative permittivity by simply assigning to the
⑀
1
and
⑀
2
coefficients the values of the relative permittivity at the
discrete points is a poor representation along the bound-
aries inside the grating region. Instead, we use an inter-
polation scheme that locally averages the permittivity or
its inverse according to the following rules:
⑀
1
共
i, i
兲
⫽
1
具
1/
⑀
共
x
兲
典
关
x
i⫺1
⬘
;x
i
⬘
兴
, (7a)
⑀
2
共
i, i
兲
⫽
具
⑀
共
x
兲
典
关
x
i
;x
i⫹1
兴
. (7b)
In Eqs. (7a) and (7b) the brackets indicate the intervals
over which the averaging has to be performed. The use
of this interpolation scheme has a drastic impact on the
convergence performance of the present method and is
justified in Appendix A. The differential operator D
1
and
D
2
are given by
D
1
⫽
冤
h
1
⫺1
共
1
兲
h
1
⫺1
共
2
兲
•
•
•
•
h
1
⫺1
共
N
兲
冥
⫻
冤
⫺11
⫺1 •
••
••
••
⫺11
␣
⫺1
冥
, (8a)
D
2
⫽
冤
h
2
⫺1
共
1
兲
h
2
⫺1
共
2
兲
•
•
•
•
h
2
⫺1
共
N
兲
冥
⫻
冤
1 ⫺
␣
⫺11
••
••
••
• 1
⫺11
冥
. (8b)
Fig. 1. Parameter definition for the classical grating diffraction
problems considered in this paper.
1034 J. Opt. Soc. Am. A / Vol. 17, No. 6/June 2000 P. Lalanne and J.-P. Hugonin
D
1
and D
2
are calculated on the crosses and circles, re-
spectively, in Fig. 2 The coefficient
␣
equal to
exp(⫺jk
0
n
1
⌳ sin
) at the lower left and upper right in ma-
trix D
1
and D
2
comes from the pseudoperiodicity condi-
tion for the electromagnetic fields in the grating region:
⌿(x ⫹ ⌳, z) ⫽ ⌿(x, z)exp(⫺jk
0
n
1
⌳ sin
). The vectors
h
1
and h
2
are defined by h
1
(i) ⫽ x
i⫹1
⫺ x
i
and h
2
(i)
⫽ x
i
⬘
⫺ x
i⫺1
⬘
.
The magnetic field ⌿ in the grating region is given by
⌿
共
z
兲
⫽
兺
m⫽1
N
W
m
兵
c
m
⫹
exp
共
⫺
m
z
兲
⫹ c
m
⫺
exp
关
m
共
z ⫺ h
兲
兴
其
,
(9)
where W
m
and
m
are the vector of the eigenvector matrix
W and the positive square root of the eigenvalues of the
matrix E, respectively. The c
m
⫹
and c
m
⫺
coefficients, like
the R
i
and T
i
coefficients, are unknowns.
All the unknown coefficients have to be evaluated by
matching the tangential field components ⌿ and (1/
⑀
)
⫻ (
⌿/
z) at the grating interfaces, z ⫽ 0 and z ⫽ h.
For this purpose, we first note that the fields in the inci-
dent medium and in the substrate are expressed in the
Rayleigh basis, whereas in the grating region, the field is
just known at some discrete locations, the mesh points
x
i
’s. A natural choice here would be to discretize the
electromagnetic field in the incident medium and that in
the substrate at the mesh points. We rather prefer (see
Section 4) to interpret the discrete ⌿
i
’s values computed
at the x
i
’s points in the grating region as the coefficients
of the magnetic field expansion in a set of rectangle func-
tions. For this purpose, we define the function
rect
i
(x), i ⫽ 1,... N,as
rect
i
共
x
兲
⫽
再
1ifx
i
⬘
⬍
x
⬍
x
i⫹1
⬘
0 otherwise
. (10)
By matching the boundary conditions, we obtain at the
top interface z ⫽ 0,
exp
共
⫺jk
0
n
1
sin
x
兲
⫹
兺
i
R
i
exp
共
⫺jk
xi
x
兲
⬅
兺
p⫽1
N
兺
m⫽1
N
w
m,p
rect
p
共
x
兲
关
c
m
⫹
⫹ c
m
⫺
exp
共
⫺
m
h
兲
兴
, (11a)
⫺jk
0
cos
n
1
exp
共
⫺jk
0
n
1
sin
x
兲
⫹
兺
i
jk
1,zi
⑀
1
R
i
exp
共
⫺jk
xi
x
兲
⬅
兺
p⫽1
N
兺
m⫽1
N
m
v
m,p
rect
p
共
x
兲
关
⫺c
m
⫹
⫹ c
m
⫺
exp
共
⫺
m
h
兲
兴
,
(11b)
and at the substrate interface z ⫽ h,
兺
i
T
i
exp
共
⫺jk
xi
x
兲
⬅
兺
p⫽1
N
兺
m⫽1
N
w
m,p
rect
p
共
x
兲
关
c
m
⫹
exp
共
⫺
m
h
兲
⫹ c
m
⫺
兴
, (12a)
兺
i
⫺jk
3,zi
⑀
3
T
i
exp
共
⫺jk
xi
x
兲
⬅
兺
p⫽1
N
兺
m⫽1
N
m
v
m,p
rect
p
共
x
兲
关
⫺c
m
⫹
exp
共
⫺
m
h
兲
⫹ c
m
⫺
兴
.
(12b)
In Eqs. (11) and (12) the sign ⬅ is used to specify the
equality between two functions of the variable x that are
expressed in different function expansions. w
m,p
is the
pth component of vector W
m
and v
m,p
is the pth compo-
nent of vector V
m
,V
m
⫽
⑀
1
⫺1
W
m
Q, with Q being a diago-
nal matrix with the diagonal element
m
. Hereafter, an
equal number of Rayleigh orders and of point locations is
considered. Thus in Eqs. (11) and (12) the summation
over the index i runs from ⫺(N ⫺ 1)/2to(N ⫺ 1)/2 for
odd N values or from ⫺N/2 to N/2 ⫺ 1 for even N values.
Equations (11) and (12) constitute a system of equations
in known function expansions with unknown expansion
coefficients R
i
, T
i
, c
m
⫹
, and c
m
⫺
. To solve this system, we
use the method of moments. In this method, a projection
basis is first chosen. Then both sides of the series-
expansion equations are projected on the projection basis.
Finally, the linear system of equations is solved with
standard numerical techniques. As the projection basis,
we chose the set of plane waves of the Rayleigh expan-
sion. This choice is commented on in Section 4. It is
easily found that rect
i
(x) ⬅ 兺
m⫽1
N
p
i,m
exp(⫺jk
xm
x), with
p
i,m
⫽ exp
共
jk
xm
x
i
兲
兵
exp
关
jk
xm
共
x
i
⬘
⫺ x
i
兲
兴
⫺ exp
关
jk
xm
共
x
i⫺1
⬘
⫺ x
i
兲
兴
其
/
共
jk
xm
⌳
兲
. (13)
In a matrix format, Eqs. (11) and (12) can be written as
冋
␦
i,0
j
␦
i,0
k
0
cos
/n
1
册
⫹
冋
I
⫺jZ
1
册
R ⫽
冋
PW PWX
PV ⫺PVX
册
冋
c
⫹
c
⫺
册
,
(14)
冋
I
jZ
2
册
T ⫽
冋
PWX PW
PVX ⫺PV
册
冋
c
⫹
c
⫺
册
, (15)
where X, Z
1
, and Z
3
are diagonal matrices with the diag-
onal elements exp(⫺
m
h), k
1,zi
/
⑀
1
and k
3,zi
/
⑀
3
, respec-
tively. P is the matrix formed by the p
i,m
coefficients.
Equations (14) and (15) are solved numerically. One
may first analytically eliminate R
i
and T
i
, then solve the
resulting set of equations for the c
m
⫹
and c
m
⫺
coefficients,
Fig. 2. Grating period discretization.
P. Lalanne and J.-P. Hugonin Vol. 17, No. 6/June 2000/J. Opt. Soc. Am. A 1035
and, finally, substitute the c
m
⫹
and c
m
⫺
coefficients back
into Eqs. (14) and (15). The diffraction efficiencies are
DE
ri
⫽ R
i
R
i
*
Re
共
k
1.zi
/k
0
n
1
cos
兲
,
DE
ti
⫽ T
i
T
i
*
Re
共
k
3,zi
n
1
/k
0
n
3
2
cos
兲
.
The field inside the grating region can be computed at the
mesh points with Eq. (9).
B. TE Polarization
We briefly outline here the few modifications necessary to
implement the TE polarization case. Equations (1)–(3)
remain valid, but now, ⌿ represent the electric field.
Discretizing the Helmholtz equation (
2
⌿)/(
z
2
)
⫹ (
2
⌿)/(
x
2
) ⫹ k
0
2
⑀
⌿ ⫽ 0 leads to
2
⌿
z
2
⫹ E⌿ ⫽ 0,
where E ⫽ D
2
D
1
⫹ k
0
2
⑀
3
. (16)
In Eq. (16),
⑀
3
is the diagonal matrix with diagonal coef-
ficients
⑀
3
共
i,i
兲
⫽
具
⑀
共
x
兲
典
关
x
i⫺1
⬘
;x
i
⬘
兴
. (17)
For TE polarization, the tangential field components are
⌿ and (
⌿)/(
z). Equations (14) and (15) become
冋
␦
i,0
j
␦
i,0
k
0
n
1
cos
册
⫹
冋
I
⫺jY
1
册
R ⫽
冋
PW PWX
PV ⫺PVX
册
冋
c
⫹
c
⫺
册
,
(18)
冋
I
jY
2
册
T ⫽
冋
PWX PW
PVX ⫺PV
册
冋
c
⫹
c
⫺
册
, (19)
where Y
1
and Y
3
are diagonal matrices with the diagonal
elements k
1,zi
and k
3,zi
, respectively. V is the matrix
formed by the vectors V
m
, V
m
⫽ W
m
Q. W, Q, and X are
as defined previously. The diffraction efficiencies are
DE
ri
⫽ R
i
R
i
*
Re
共
k
1,zi
/k
0
n
1
cos
兲
,
DE
ti
⫽ T
i
T
i
*
Re
共
k
3,zi
/k
0
n
1
cos
兲
.
3. NUMERICAL RESULTS
In this section the convergence performance of the
present method is studied for TM and TE polarizations.
For the sake of comparison, the convergence speed
achieved is systematically compared with that achieved
with the RCWA.
26
Our RCWA implementation for TM
polarization is described in Refs. 7 and 8. As will be
shown, the relative performance with respect to the
RCWA strongly depends on the grating materials. Me-
tallic lamellar gratings operating in the visible and in the
infrared region of the spectrum are first considered in
Subsection 3.A. Dielectric gratings for which relatively
bad convergence performance are observed are then con-
sidered in Subsection 3.B. In all calculations, the num-
ber of retained Rayleigh orders is equal to the number N
of point locations.
A. Metallic Gratings
In general, with FD methods, the convergence perfor-
mance depends on the set of selected discrete points. We
first consider a uniform sampling and a metallic grating
on an aluminum substrate in a Littrow mounting for op-
eration with CO
2
lasers. This example was previously
considered in the literature by several researchers
7,8,23
concerned with the convergence performance of the
RCWA. The diffraction parameters are given in the cap-
tion of Fig. 3. Only the minus-first and zeroth reflected
orders are propagating. For a perfectly uniform discreti-
zation, Figs. 3(a) and 3(b) show the error (defined as the
computed diffraction efficiency minus the exact value) of
the zeroth order as a function of N for TM and TE polar-
izations, respectively. In these figures, pluses are ob-
tained with the present method and circles with the
RCWA. This convention is used throughout the paper.
Fig. 3. Uniform discretization: computational error for the re-
flected zeroth order of a metallic grating: (a) TM, (b) TE. The
values of the parameters are:
⫽ 30°, f ⫽ 0.5, n
1
⫽ n
⬘
⫽ 1,
n
3
⫽ n ⫽ 0.22 ⫹ 6.71i, and ⫽ ⌳ ⫽ h ⫽ 1.0
m. The exact
value of the reflected zeroth order is 84.848% for TM and
73.428% for TE. Pluses, present method; circles, RCWA.
1036 J. Opt. Soc. Am. A / Vol. 17, No. 6/June 2000 P. Lalanne and J.-P. Hugonin
