Adaptive perfectly matched layer for Wood’s anomalies in diffraction gratings

doi:10.1364/OE.20.028094

« Show journal navigation

Adaptive perfectly matched layer for Wood’s anomalies in diffraction gratings

Benjamin Vial, Frédéric Zolla, André Nicolet, Mireille Commandré, and Stéphane Tisserand »View Author Affiliations

Optics Express, Vol. 20, Issue 27, pp. 28094-28105 (2012)
http://dx.doi.org/10.1364/OE.20.028094

View Full Text Article

Acrobat PDF (1269 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Setup of the problem and notations
3. Construction of an adaptive PML
4. Numerical example
5. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (16)
Cited By
Figures (7)
Metrics

Abstract

We propose an Adaptive Perfectly Matched Layer (APML) to be used in diffraction grating modeling. With a properly tailored co-ordinate stretching depending both on the incident field and on grating parameters, the APML may efficiently absorb diffracted orders near grazing angles (the so-called Wood’s anomalies). The new design is implemented in a finite element method (FEM) scheme and applied on a numerical example of a dielectric slit grating. Its performances are compared with classical PML with constant stretching coefficient.

1. Introduction

Since their introduction by Bérenger in [1

J.-P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics 114, 185 – 200 (1994). [CrossRef]

] for the time dependent Maxwell’s equations, Perfectly Matched Layers (PMLs) have become a widely used technique in computational physics. The idea is to enclose the area of interest by surrounded layers which are absorbing and perfectly reflectionless. Subsequent formulations [2

M. Lassas, J. Liukkonen, and E. Somersalo, “Complex riemannian metric and absorbing boundary conditions,” J. Math. Pure Appl. 80, 739 – 768 (2001).

, 3

M. Lassas and E. Somersalo, “Analysis of the PML equations in general convex geometry,” P. Roy. Soc. Edinb. A 131, 1183–1207 (2001). [CrossRef]

] allow to consider PMLs as a complex change of co-ordinates, which implies equivalent material properties [4

A. Nicolet, F. Zolla, Y. Ould Agha, and S. Guenneau, “Geometrical transformations and equivalent materials in computational electromagnetism,” Compel 27, 806–819 (2008). [CrossRef]

, 5

F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau, D. Felbacq, A. Argyros, and S. Leon-Saval, Foundations of Photonic Crystal Fibres (Imperial College Press, London, 2012), 2nd ed.

]. A properly designed change of variables preserves the solution in the region of interest while introducing exponential decay at infinity. It is then obvious to truncate the problem to a finite domain, set a convenient boundary condition on the boundary of this domain and apply the FEM. In a scattering problem, the “traditional co-ordinate stretching” is frequency dependent, yielding a constant attenuation of propagating waves. However, these kinds of PMLs are inefficient for periodic problems when dealing with grazing angles of diffracted orders, i.e. when the frequency is near a Wood’s anomaly [6

R. Wood, “On a remarkable case of uneven distribution of light in a diffraction grating spectrum,” P. Phys. Soc. Lond. 18, 269–275 (1902). [CrossRef]

, 7

Lord Rayleigh, “Note on the remarkable case of diffraction spectra described by prof. Wood,” Philos. Mag. 14, 60–65 (1907). [CrossRef]

], leading to spurious reflexions and thus numerical pollution of the results. An important question in designing absorbing layers is thus the choice of their parameters : the PML thickness and the absorption coefficient. To this aim, adaptive formulation have been set up, most of them employing a posteriori error estimate [8

Z. Chen and X. Liu, “An adaptive perfectly matched layer technique for time-harmonic scattering problems,” SIAM J. Numer. Anal. 43, 645–671 (2005). [CrossRef]

, 9

G. Bao, Z. Chen, and H. Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A 22, 1106–1114 (2005). [CrossRef]

, 10

A. Schädle, L. Zschiedrich, S. Burger, R. Klose, and F. Schmidt, “Domain decomposition method for Maxwell’s equations: scattering off periodic structures,” J. Comput. Phys. 226, 477 – 493 (2007). [CrossRef]

]. In this paper, we propose Adaptive PMLs (APMLs) with a suitable co-ordinate stretching, depending both on incidence and grating parameters, capable of efficiently absorbing propagating waves with nearly grazing angles. In the first section, we detail our FEM formulation and expose the problems arising when dealing with Wood’s anomalies. Section 2 is dedicated to the mathematical formulation used to determine PML parameters adapted to any diffraction orders. In the last section, we provide a numerical example of a dielectric slit grating showing the efficiency of our approach in comparison with classical PMLs.

2. Setup of the problem and notations

2.1. Governing equations and description of the structure

We denote by x, y and z the unit vectors of an orthogonal co-ordinate system Oxyz. We deal with time-harmonic fields, so that the electric and magnetic fields are represented by complex vector fields E and H with a time-dependence in exp(−iωt), which will be discarded in the sequel. Moreover, we denote k₀ = ω/c.

The materials in this paper may be z-anisotropic, so the tensor fields of relative permittivity ε͇ and relative permeability μ͇ are of the following form :

\underset{=}{ε} = (\begin{matrix} ε_{x x} & \bar{ε_{a}} & 0 \\ ε_{a} & ε_{y y} & 0 \\ 0 & 0 & ε_{z z} \end{matrix}) and \underset{=}{μ} = (\begin{matrix} μ_{x x} & \bar{μ_{a}} & 0 \\ μ_{a} & μ_{y y} & 0 \\ 0 & 0 & μ_{z z} \end{matrix}),

(1)

where the coefficients ε_xx, ε_aa,...μ_zz are (possibly) complex valued functions of x and y, and where

\bar{ε_{a}}

(resp.

\bar{μ_{a}}

) is the complex conjugate of ε_a (resp. μ_a).

The grating is illuminated by an incident plane wave of wave vector defined by the angle θ₀ : k⁺ = α⁺x+β⁺y = k⁺(sinθ₀x − cosθ₀y). Its electric (resp. magnetic) field is linearly polarized along the z-axis, this is the so-called transverse electric or s-polarization case (resp. transverse magnetic or p-polarization case) :

E^{0} = A_{e}^{0} e^{(i k^{+} \cdot r)} z (resp . H^{0} = A_{m}^{0} e^{(i k^{+} \cdot r)} z)

(2)

where

A_{e}^{0}

(resp.

A_{m}^{0}

) are arbitrary complex numbers and r = (x,y)^T.

The diffraction problem we are dealing with is to solve Maxwell’s equations in harmonic regime, i.e. to find the unique solution (E, H) of :

{\begin{cases} curl E = i ω μ_{0} \underset{=}{μ} H \\ curl H = - i ω ε_{0} \underset{=}{ε} E, \end{cases}

(3a)(3b)

such that the diffracted field E^d := E − E⁰ and H^d := H − H⁰ satisfies an Outgoing Wave Condition (OWC, [11

A. Sommerfeld, Partial Differential Equations in Physics (Academic Press, New York, 1949).

]) and where E and H are quasi-periodic with respect to the x co-ordinate.

Under the aforementioned assumptions, the diffraction problem in non conical mounting can be separated in two fundamental scalar cases TE and TM. Thus we search for a z-linearly polarized electric (resp.magnetic) field E = e(x,y)z (resp. H = h(x,y)z). Denoting ε͇̃ and μ͇̃ the 2 × 2 matrices extracted from ε͇ and μ͇:

\underset{=}{\tilde{ε}} = (\begin{matrix} ε_{x x} & \bar{ε_{a}} \\ ε_{a} & ε_{y y} \end{matrix}) and \underset{=}{\tilde{μ}} = (\begin{matrix} μ_{x x} & \bar{\bar{μ_{a}}} \\ μ_{a} & μ_{y y} \end{matrix}),

(4)

the functions e and h are solution of similar differential equations :

ℒ_{\underset{=}{ξ}, χ} (u) : = div (\underset{=}{ξ} \cdot grad u) + k_{0}^{2} χ u = 0,

(5)

with

u = e, \underset{=}{ξ} = {\underset{=}{\tilde{μ}}}^{T} / det (\underset{=}{\tilde{μ}}), χ = ε_{z z},

(6)

for the TE case,

u = h, \underset{=}{ξ} = {\underset{=}{\tilde{ε}}}^{T} / det (\underset{=}{\tilde{ε}}), χ = μ_{z z},

(7)

for the TM case.

The gratings we are dealing with are made of three regions (see Fig. 1) :

The superstrate (y > h^g) which is supposed to be homogeneous, isotropic and lossless and characterized solely by its relative permittivity ε⁺ and its relative permeability μ⁺ and we denote $k^{+} : = k_{0} \sqrt{ε^{+} μ^{+}}$ , where k₀ := ω/c.
The substrate (y < 0) which is supposed to be homogeneous and isotropic and therefore characterized by its relative permittivity ε⁻ and its relative permeability μ⁻ and we denote $k^{-} : = k_{0} \sqrt{ε^{-} μ^{-}}$
The groove region (0 < y < h^g), embedded in the superstrate, which can be heterogeneous, z-anisotropic and lossy, so that it is characterised by the tensor fields ε͇_g(x,y) and μ͇_g(x,y). The periodicity of the grooves along Ox is denoted d.

Fig. 1 Set up of the problem and notations

2.2. Implementation of the PMLs

Transformation optics have recently unified various techniques in computational electromagnetics such as the treatment of open problems, helicoidal geometries or the design of invisibility cloaks [4

A. Nicolet, F. Zolla, Y. Ould Agha, and S. Guenneau, “Geometrical transformations and equivalent materials in computational electromagnetism,” Compel 27, 806–819 (2008). [CrossRef]

]. These apparently different problems share the same concept of geometrical transformation, leading to equivalent material properties. A very simple and practical rule can be set up [5

F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau, D. Felbacq, A. Argyros, and S. Leon-Saval, Foundations of Photonic Crystal Fibres (Imperial College Press, London, 2012), 2nd ed.

] : when changing the co-ordinate system, all you have to do is to replace the initial materials properties ε͇ and μ͇ by equivalent material properties ε͇_s and μ͇_s given by the following rule :

{\underset{=}{ε}}_{s} = J^{- 1} \underset{=}{ε} J^{- T} det (J) and {\underset{=}{μ}}_{s} = J^{- 1} \underset{=}{μ} J^{- T} det (J),

(8)

where J is the Jacobian matrix of the co-ordinate transformation consisting of the partial derivatives of the new co-ordinates with respect to the original ones (J^−T is the transposed of its inverse).

In this framework, the most natural way to define PMLs is to consider them as maps on a complex space ℂ³, which co-ordinate change leads to equivalent permittivity and permeability tensors. We detail here the different co-ordinates used in this paper.

(x,y,z) are the cartesian original co-ordinates.
(x_s, y_s, z_s) are the complex stretched co-ordinates. A suitable subspace Γ ⊂ ℂ³ is chosen (with three real dimensions) such that (x_s, y_s, z_s) are the complex valued co-ordinates of a point on Γ (e.g. x = 𝔕𝔢(x_s), y = 𝔕𝔢(y_s), z = 𝔕𝔢(z_s)).
(x_c, y_c, z_c) are three real co-ordinates corresponding to a real valued parametrization of Γ ⊂ ℂ³.

In this paper, we use rectangular PMLs [12

Y. Ould Agha, A. Nicolet, F. Zolla, and S. Guenneau, “Leaky modes in twisted microstructured optical fibres,” Wave Random Complex 17, 559–570 (2007). [CrossRef]

] absorbing in the y-direction and we choose a diagonal matrix J = diag(1, s_y(y), 1), where s_y(y) is a complex-valued function defined by :

y_{s} (y) = \int_{0}^{y} s_{y} (y^{'}) d y^{'} .

(9)

The expression of the equivalent permittivity and permeability tensors are thus :

{\underset{=}{ε}}_{s} = (\begin{matrix} s_{y} ε_{x x} & \bar{ε_{a}} & 0 \\ ε_{a} & s_{y}^{- 1} ε_{y y} & 0 \\ 0 & 0 & s_{y} ε_{z z} \end{matrix}) and {\underset{=}{μ}}_{s} = (\begin{matrix} s_{y} μ_{x x} & \bar{μ_{a}} & 0 \\ μ_{a} & s_{y}^{- 1} μ_{y y} & 0 \\ 0 & 0 & s_{y} μ_{z z} \end{matrix}) .

(10)

Note that the equivalent medium has the same impedance than the original one as ε͇ an μ͇ are transformed in the same way, which guarantees that the PML is perfectly reflectionless.

We are now in position to define the so-called substituted field F_s = (E_s, H_s), solution of Eqs. (3) with ξ͇ = ξ͇_s and χ = χ_s. F_s equals the field F in the region y^b < y < y^t (with y^b = −h⁻ and y^t = h^g + h⁺), provided that s_y(y) = 1 in this region (cf. Fig. 2). The complex co-ordinate mapping y(y_c) is simply defined as the derivative of the stretching coefficient s_y(y) with respect to y_c. With simple stretching functions, we can obtain a reliable criterion on the decaying of the fields. A classical choice is :

s_{y} (y) = {\begin{array}{l} ζ^{-} & if y < y^{b} \\ 1 & if y^{b} < y < y^{t} \\ ζ^{+} & if y > y^{t} \end{array}

(11)

where ζ^± = ζ^′,± + iζ″^,± are complex constants with ζ″^,± > 0.

Fig. 2 The basic cell used for the FEM computation of the diffracted field

u_{2}^{d}

In that case, the complex valued function y(y_c) defined by Eq. (9) is explicitly given by :

y (y_{c}) = {\begin{array}{l} y^{b} + ζ^{-} (y_{c} - y^{b}) & if y_{c} < y^{b} \\ y_{c} & if y^{b} < y_{c} < y^{t}, \\ y^{t} + ζ^{+} (y_{c} - y^{t}) & if y_{c} > y^{t} \end{array}

(12)

Let’s now study the behaviour of the field in the PMLs. In the substrate and the superstrate, according to Bloch’s theorem, the diffracted field u^d(x,y) can be written as :

u^{d} (x, y) = \sum_{n \in ℤ} u_{n}^{d} (y) e^{i α_{n} x},

(13)

where

u_{n}^{d} (y) = \frac{1}{d} \int_{- d / 2}^{d / 2} u^{d} (x, y) e^{- i α_{n} x} d x, with α_{n} = α + \frac{2 π}{d} n .

(14)

Eventually, denoting

β_{n}^{\pm^{2}} = k^{\pm^{2}} - α_{n}^{\pm^{2}}

(15)

and accounting for an OWC, we obtain the expressions of the different diffraction orders :

u_{n}^{d} (y) = {\begin{matrix} u_{n}^{+} (y) = r_{n} e^{i β_{n}^{+} y} & for y > h^{g} \\ u_{n}^{-} (y) = t_{n} e^{- i β_{n}^{-} y} & for y < 0 . \end{matrix}

(16)

We focus on the bottom PML, as similar considerations can be conducted for the top PML. We consider a diffraction order propagating in the substrate, which can be expressed in the stretched co-ordinate system as :

u_{n, s}^{-} (y_{c}) = u_{n}^{-} (y (y_{c})) = t_{n} e^{- i β_{n}^{-} [y^{t} + ζ^{-} (y_{c} - y^{t})]} .

The non oscillating part of this function is given by :

U_{n}^{-} (y) = t_{n} exp ((β_{n}^{', -} ζ^{″, -} + β_{n}^{″, -} ζ^{', -}) y_{c}),

where

β_{n}^{-} = β_{n}^{', -} + i β_{n}^{″, -}

. For a propagating order we have

β_{n}^{', -} > 0

and

β_{n}^{″, -} = 0

, while for an evanescent order

β_{n}^{', -} = 0

and

β_{n}^{″, -} > 0

. It is thus sufficient to take ζ^′,− > 0 and ζ″^,− > 0 to ensure the exponential decay to zero of the field inside the PML if it was of infinite extent. But, of course, for practical purposes, the height of the PML is finite and have to be suitably chosen. For this, two pitfalls must be avoided :

The PML is too small compared to the skin depth. As a consequence, the electromagnetic wave cannot be considered as vanishing : this limited PML is no longer reflectionless.
The PML is much larger than the skin depth. In that case, a significant part of the PML is not useful, which gives rise to the resolution of linear systems of large dimensions.

It then remains to derive the skin depth,

l_{n}^{-}

, associated with the propagating order n. This characteristic depth is defined as the length for which the field is divided by e :

U_{n}^{-} (y - l_{n}^{-}) = \frac{U_{n}^{-} (y)}{e} .

l_{n}^{-} = {(β_{n}^{', -} ζ^{″, -} + β_{n}^{″, -} ζ^{', -})}^{- 1}

, and we take the larger value among the

l_{n}^{-}

l^{-} = max_{n \in ℤ} l_{n}^{-} .

We set the height of the bottom PML region ĥ⁻ = 10l⁻.

2.3. The FEM formulation

Our FEM formulation is that described in [13

G. Demésy, F. Zolla, A. Nicolet, M. Commandré, and C. Fossati, “The finite element method as applied to the diffraction by an anisotropic grating,” Opt. Express 15, 18089–18102 (2007). [CrossRef] [PubMed]

, 14

G. Demésy, F. Zolla, A. Nicolet, and M. Commandré, “Versatile full-vectorial finite element model for crossed gratings,” Opt. Lett. 34, 2216–2218 (2009). [CrossRef] [PubMed]

]. It relies on the fact that the diffraction problem can be rigorously treated as an equivalent radiation problem with sources inside the diffractive object. The expression of the source terms depends on the field u₁ solution of the annex problem of a simple interface, and is known in closed form. In this context, the role of the PMLs is to damp the field radiated by the sources. We set quasi-periodicity conditions on lateral boundaries and Neumann homogeneous conditions are imposed on the outward boundary of each PML. By so doing, the electromagnetic wave is not forced to vanish and the values of the computed field on these boundaries give valuable information on the absorbing efficiency of the PMLs. The cell defined in Fig. 2 is meshed using 2^nd order Lagrange elements [15

A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements (Springer, New York, 2004).

]. In the numerical examples in the sequel, the maximum element size is set to

λ_{0} / (N_{m} \sqrt{ℜ 𝔢 (ε_{r})})

, where N_m is an integer (between 6 and 10 is usually a good choice). The final algebraic system is solved using a direct solver (PARDISO).

Note that we use this FEM formulation in this paper but the results of the sequel about PMLs are general and can be applied to any numerical method used for the modelling of diffraction gratings.

2.4. Weakness of the classical PML for grazing diffracted angles

Resuming part 2.2, we consider only the bottom PML, as analogous conclusions holds for the top PML. The efficiency of the classical PML fails for grazing diffracted angles, in other words when a given order appears/vanishes : this is the so-called Wood’s anomaly, well known in the theory of gratings. In mathematical words, there exists n₀ such that

β_{n_{0}}^{-} ≃ 0

. The skin depth of the PML then becomes very large. To compensate it, one should increase the value of ζ^″,−, but this gives rise to spurious numerical reflections due to an overdamping. For a fixed value of ĥ⁻, if ζ^″^,− is too weak, the absorption in the PMLs is insufficient and the wave reflects on the outward boundary of the PML. To illustrate these typical behaviours (cf. Fig. 3), we compute the diffracted field by a grating with a rectangular cross section of height h^g = 1.5μm and width L^g = 3μm with ε^g = 11.7, deposited on a substrate with permittivity ε⁻ = 2.25. The structure is illuminated by a p-polarized plane wave of wavelength λ₀ = 10μm and of angle of incidence θ₀ = 10° in the air (ε⁺ = 1). All materials are non magnetic (μ_r = 1) and the periodicity of the grating is d = 4μm. We set

{\hat{h}}^{-} = 10 l_{0}^{-}

and ζ^′,− = 1.

Fig. 3 Zero^th transmitted order by a grating with a rectangular cross section (see parameters in text, part 2.4) for different values of ζ^″,− : blue line, ζ^″,− = 1, correct damping; green line, ζ^″,− = 0.1, underdamping; red line, ζ^″,− = 20, overdamping.

3. Construction of an adaptive PML

To overcome the problems pointed out in the preceding section, we propose a co-ordinate stretching that rigorously treats the problem of Wood’s anomalies. The wavelengths “seen” by the system are very different depending on the order at stake :

if the diffracted angle θ_n is zero, the apparent wavelength λ/cosθ_n is simply the incident wavelength,
if the diffracted angle is near ±π/2 (grazing angle), the apparent wavelength λ/cosθ_n is very large.

Thus if a classical PML is adapted to one diffracted order, it will not be for another, and vice versa. The idea behind the APML is to deal with each and every order when progressing in the absorbing medium.

Once again the development will be conducted only for the PML adapted to the substrate. We consider a real-valued co-ordinate mapping y_d(y), the final complex-valued mapping is then y_c(y) = ζ⁻y_d(y), with the complex constant ζ⁻, with ζ^′,− > 0 and ζ^″,− > 0, accounting for the damping of the PML medium. We start by transforming Eq. (15),

n \mapsto β_{n}^{-}

a function with integer argument into a function with real argument continuously interpolated between the imposed integer values. Indeed, the geometric transformations associated to the PML has to be continuous and differentiable in order to compute the Jacobian. For this purpose, we choose the parametrization :

α (y_{d}) = α_{0} + \frac{2 π}{d} \frac{y_{d}}{λ_{0}},

(17)

so that the application β⁻ defined by

β^{-} {(y_{d})}^{2} = k_{0}^{2} ε^{-} - α {(y_{d})}^{2}

is continuous. Thus, the propagation constant of the n^th transmitted order is given by

β_{n}^{-} = β^{-} (n λ_{0})

. The key idea is to combine the complex stretching with a real non uniform contraction (given by the continuous function y(y_d), Eq.(19)). This contraction is chosen in such a way that for each order n there is a depth

y_{d}^{n}

such that around this depth the apparent wavelength corresponding to the order in play is contracted to a value close to λ₀. At that point of the PML, this order is perfectly absorbed thanks to the complex stretch. We thus eliminate first the orders with quasi normal diffracted angles at lowest depths up to grazing orders (near Wood’s anomalies) which are absorbed at greater depths. In mathematical words, the translation of previous considerations on the real contraction can be expressed as :

exp [- i β^{-} (y_{d}) y (y_{d})] = exp (- i k_{0} y_{d})

(18)

The contraction y(y_d) is thus given by :

y (y_{d}) = \frac{k_{0} y_{d}}{β^{-} (y_{d})} = \frac{y_{d}}{\sqrt{ε^{-} - {(sin θ_{0} + y_{d} / d)}^{2}}}

(19)

The function y(y_d) has two poles, denoted

y_{d, \pm}^{★} = d (\pm \sqrt{ε^{-}} - sin θ_{0})

. When

y_{d, \pm}^{★} = \pm n λ_{0}

with n ∈ ℕ^★,

β^{-} (y_{d, \pm}^{★}) = β^{-} (\pm n λ_{0}) = β_{\pm}^{-} = 0

, i.e. we are on a Wood’s anomaly associated with the appearance/disappearance of the ±n^th transmitted order. We now search for the nearest point to

y_{d, \pm}^{*}

associated with a Wood’s anomaly, denoting :

{\begin{matrix} n_{+}^{★} / & D_{+} = min_{n_{+}^{★} \in ℕ^{★}} | y_{d, +}^{★} - n_{+}^{★} λ_{0} | \\ n_{-}^{★} / & D_{-} = min_{n_{-}^{★} \in ℕ^{★}} | y_{d, -}^{★} + n_{-}^{★} λ_{0} |, \end{matrix}

In a second step, we seek the point

y_{d}^{0} = n^{★} λ_{0}

such that :

n^{★} / D = min_{n^{★} \in {n_{+}^{★}, n_{-}^{★}}} (D^{+}, D^{-}),

(20)

To avoid the singular behaviour at

y_{d} = y_{d, \pm}^{★}

, we continue the graph of the function y_d(y) by a straight line tangent at

y_{d}^{0}

, which equation is

t_{0} (y_{d}) = s (y_{d}^{0}) (y_{d} - y_{d}^{0}) + y (y_{d}^{0})

, where

s (y_{d}) = \frac{\partial y}{\partial y_{d}} (y_{d})

is the so-called stretching coefficient. The final change of co-ordinate is then given by :

\tilde{y} (y_{d}) = {\begin{matrix} y (y_{d}) & for y_{d} \leq y_{d}^{0} \\ t_{0} (y_{d}) & for y_{d} > y_{d}^{0} . \end{matrix}

(21)

Figure 4 shows an example of this co-ordinate mapping.

Fig. 4 Example of a co-ordinate mapping ỹ(y_d) used for the APML (black solid line). The graph of y_d(y) (blue solid line) is continued by a straight line t₀(y_d) tangent at

y_{d}^{0}

(red dashed line) to avoid the singular behaviour at

y_{d} = y_{d}^{★}

Eventually, the complex stretch s_y used in Eq. (10) is given by :

s_{y} (y_{d}) = ζ^{-} \frac{\partial \tilde{y}}{\partial y_{d}} (y_{d}) .

(22)

Equipped with this mathematical formulation, we can tailor a layer that is doubly perfectly matched :

to a given medium, which is the aim of the PML technique, through Eq. (8),
to all diffraction orders, through the stretching coefficient s_y, which depends on the characteristics of the incident wave and on opto-geometric parameters of the grating.

4. Numerical example

We now apply the method described in the preceding parts to design an adapted bottom PML for the same example as in part 2.4. The parameters are the same, and we choose the wavelength of the incident plane wave close to the Wood’s anomaly related to the +1 transmitted order (

λ_{0} = 0.999 y_{d, +}^{★}

). Moreover, we set the length of the PML

{\hat{h}}^{-} = 1.1 y_{d, +}^{★}

and choose absorption coefficients ζ⁺ = ζ⁻ = 1 + i. For both cases (PML and APML), parameters are alike, the only difference being the complex stretch s_y.

The field maps of the norm of H_z, E_x and E_y are plotted in logarithmic scale on Fig. 5, for the case of a classical PML and our APML. We can observe that the field H_z that is effectively computed is clearly damped in the bottom APML (leftmost on Fig. 5(b)) whereas it is not in the standard case (leftmost on Fig. 5(a)), causing spurious reflections on the outer boundary. The fields E_x and E_y are deduced from H_z thanks to Maxwell’s equations. The high values of E_y at the tip of the APML (rightmost on Fig. 5(b)) are due to very high values of the optical equivalent properties of the APML medium (due to high values of s_y), which does not affect the accuracy of the computed field within the domain of interest.

Fig. 5 Field maps of the logarithm of the norm of H_z, E_x and E_y for the dielectric slit grating at

λ_{0} = 0.999 y_{d, +}^{★}

(same parameters as in part 2.4). (a) : classical PML with inefficient damping of H_z in the bottom PML. (b) : APML where the H_z field is correctly damped in the bottom PML. For both cases the thickness of the PML is

{\hat{h}}^{-} = 1.1 y_{d, +}^{★}

Another feature of our approach is that it efficiently absorbs the grazing diffraction order, as illustrated on Fig. 6 : the +1 transmitted order does not decrease in the standard PML (blue solid line), and reaches a high value at y = −ĥ⁻, whereas the same order tends to zero as y → −ĥ⁻ in the case of the adapted PML (blue dashed line).

Fig. 6 Modulus of the u_n for the three propagating orders with adapted (dashed lines) and classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except for the grazing one (n = 1) as expected. This drawback is bypassed when using the adaptive PML.

To further validate the accuracy of the method, we compare the diffraction efficiencies computed by our FEM formulation with PML and APML to those obtained by another method. We choose the Rigorous Coupled Wave Analysis (RCWA), also known as the Fourier Modal Method (FMM, [16

L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758–2767, (1997). [CrossRef]

]). For the chosen parameters, only the 0^th order is propagative in reflexion and the orders −1, 0 and +1 are non evanescent in transmission. We can also check the energy balance B = R₀ + T₋₁ + T₀ + T₊₁ since there is no lossy medium in our example. Results are reported in Table 1, and show a good agreement of the FEM with APML with the results from RCWA. On the contrary, if classical PML are used, the diffraction efficiencies are less accurate compared to those computed with RCWA. Checking the energy balance leads the same conclusions : the numerical result is perturbed by the reflection of the waves at the end of the PML if it is not adapted to the situation of nearly grazing diffracted orders.

Table 1 Diffraction efficiencies R₀, T₋₁, T₀ and T₊₁ of the four propagating orders, and energy balance B = R₀ + T₋₁ + T₀ + T₊₁, computed by three methods : RCWA (line 1), FEM formulation with APML (line 2), FEM formulation with classical PML (line 3).

	R₀	T₋₁	T₀	T₊₁	B

RCWA	0.1570	0.3966	0.1783	0.2680	0.9999
FEM + APML	0.1561	0.3959	0.1776	0.2703	0.9999
FEM + PML	0.1904	0.4118	0.1927	0.2481	1.0430

Eventually, to illustrate the behaviour of the adaptive PML when the incident wavelength gets closer to a given Wood’s anomaly, we computed the mean value of the norm of H_z along the outer boundary of the bottom PML γ = 〈|H_z(−ĥ⁻)|〉_x, when

λ_{0} = (1 + 10^{- n}) y_{d, +}^{★}

and

λ_{0} = (1 - 10^{- n}) y_{d, +}^{★}

, for n = 1, 2,...10. The results are shown in Fig. 7. As the wavelength gets closer to

y_{d, +}^{★}

, γ first increases but for n > 3, it decreases exponentially. However, in all cases, the value of γ remains small enough to ensure the efficiency of the PMLs.

Fig. 7 Mean value of the norm of H_z along the outer boundary of the bottom PML γ = 〈|H_z(−ĥ⁻)|〉_x, for λ₀ approaching the Wood’s anomaly

y_{d, +}^{★}

by inferior values (

λ_{0} = (1 - 10^{- n}) y_{d, +}^{★}

, red squares) and by superior value (

λ_{0} = (1 + 10^{- n}) y_{d, +}^{★}

, blue circles) as a function of n.

5. Conclusion

In this paper, we have proposed an adaptive PML that can treat rigorously Wood’s anomalies in numerical analysis of diffraction gratings. It is based on a complex-valued co-ordinate stretching that deals with grazing diffracted orders, yielding an efficient absorption of the field inside the PML. We provided an example in the TM polarization case (but similar results hold for the TE case), illustrating the efficiency of our method. The value of the magnetic field on the outward boundary of the PML remains small enough to consider there is no spurious reflection. The formulation is used with the FEM but can be applied to others numerical methods. Moreover, the generalization to the vectorial three-dimensional case is straightforward : the recipes given in part 3 do work irrespective of the dimension and whether the problem is vectorial. In addition, although the designed co-ordinate stretching is specific to the context of Wood’s anomalies, one can adapt this kind of non uniform PML to others wave problems by following the methodology exposed here.

References and links

1.	J.-P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics 114, 185 – 200 (1994). [CrossRef]
2.	M. Lassas, J. Liukkonen, and E. Somersalo, “Complex riemannian metric and absorbing boundary conditions,” J. Math. Pure Appl. 80, 739 – 768 (2001).
3.	M. Lassas and E. Somersalo, “Analysis of the PML equations in general convex geometry,” P. Roy. Soc. Edinb. A 131, 1183–1207 (2001). [CrossRef]
4.	A. Nicolet, F. Zolla, Y. Ould Agha, and S. Guenneau, “Geometrical transformations and equivalent materials in computational electromagnetism,” Compel 27, 806–819 (2008). [CrossRef]
5.	F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau, D. Felbacq, A. Argyros, and S. Leon-Saval, Foundations of Photonic Crystal Fibres (Imperial College Press, London, 2012), 2nd ed.
6.	R. Wood, “On a remarkable case of uneven distribution of light in a diffraction grating spectrum,” P. Phys. Soc. Lond. 18, 269–275 (1902). [CrossRef]
7.	Lord Rayleigh, “Note on the remarkable case of diffraction spectra described by prof. Wood,” Philos. Mag. 14, 60–65 (1907). [CrossRef]
8.	Z. Chen and X. Liu, “An adaptive perfectly matched layer technique for time-harmonic scattering problems,” SIAM J. Numer. Anal. 43, 645–671 (2005). [CrossRef]
9.	G. Bao, Z. Chen, and H. Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A 22, 1106–1114 (2005). [CrossRef]
10.	A. Schädle, L. Zschiedrich, S. Burger, R. Klose, and F. Schmidt, “Domain decomposition method for Maxwell’s equations: scattering off periodic structures,” J. Comput. Phys. 226, 477 – 493 (2007). [CrossRef]
11.	A. Sommerfeld, Partial Differential Equations in Physics (Academic Press, New York, 1949).
12.	Y. Ould Agha, A. Nicolet, F. Zolla, and S. Guenneau, “Leaky modes in twisted microstructured optical fibres,” Wave Random Complex 17, 559–570 (2007). [CrossRef]
13.	G. Demésy, F. Zolla, A. Nicolet, M. Commandré, and C. Fossati, “The finite element method as applied to the diffraction by an anisotropic grating,” Opt. Express 15, 18089–18102 (2007). [CrossRef] [PubMed]
14.	G. Demésy, F. Zolla, A. Nicolet, and M. Commandré, “Versatile full-vectorial finite element model for crossed gratings,” Opt. Lett. 34, 2216–2218 (2009). [CrossRef] [PubMed]
15.	A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements (Springer, New York, 2004).
16.	L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758–2767, (1997). [CrossRef]

OCIS Codes
(050.0050) Diffraction and gratings : Diffraction and gratings
(050.1755) Diffraction and gratings : Computational electromagnetic methods

ToC Category:
Diffraction and Gratings

History
Original Manuscript: September 19, 2012
Revised Manuscript: November 5, 2012
Manuscript Accepted: November 11, 2012
Published: December 4, 2012

Citation
Benjamin Vial, Frédéric Zolla, André Nicolet, Mireille Commandré, and Stéphane Tisserand, "Adaptive perfectly matched layer for Wood’s anomalies in diffraction gratings," Opt. Express 20, 28094-28105 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-27-28094

Sort: Author | Year | Journal | Reset

References

J.-P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics114, 185 – 200 (1994). [CrossRef]
M. Lassas, J. Liukkonen, and E. Somersalo, “Complex riemannian metric and absorbing boundary conditions,” J. Math. Pure Appl.80, 739 – 768 (2001).
M. Lassas and E. Somersalo, “Analysis of the PML equations in general convex geometry,” P. Roy. Soc. Edinb. A131, 1183–1207 (2001). [CrossRef]
A. Nicolet, F. Zolla, Y. Ould Agha, and S. Guenneau, “Geometrical transformations and equivalent materials in computational electromagnetism,” Compel27, 806–819 (2008). [CrossRef]
F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau, D. Felbacq, A. Argyros, and S. Leon-Saval, Foundations of Photonic Crystal Fibres (Imperial College Press, London, 2012), 2nd ed.
R. Wood, “On a remarkable case of uneven distribution of light in a diffraction grating spectrum,” P. Phys. Soc. Lond.18, 269–275 (1902). [CrossRef]
Lord Rayleigh, “Note on the remarkable case of diffraction spectra described by prof. Wood,” Philos. Mag.14, 60–65 (1907). [CrossRef]
Z. Chen and X. Liu, “An adaptive perfectly matched layer technique for time-harmonic scattering problems,” SIAM J. Numer. Anal.43, 645–671 (2005). [CrossRef]
G. Bao, Z. Chen, and H. Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A22, 1106–1114 (2005). [CrossRef]
A. Schädle, L. Zschiedrich, S. Burger, R. Klose, and F. Schmidt, “Domain decomposition method for Maxwell’s equations: scattering off periodic structures,” J. Comput. Phys.226, 477 – 493 (2007). [CrossRef]
A. Sommerfeld, Partial Differential Equations in Physics (Academic Press, New York, 1949).
Y. Ould Agha, A. Nicolet, F. Zolla, and S. Guenneau, “Leaky modes in twisted microstructured optical fibres,” Wave Random Complex17, 559–570 (2007). [CrossRef]
G. Demésy, F. Zolla, A. Nicolet, M. Commandré, and C. Fossati, “The finite element method as applied to the diffraction by an anisotropic grating,” Opt. Express15, 18089–18102 (2007). [CrossRef] [PubMed]
G. Demésy, F. Zolla, A. Nicolet, and M. Commandré, “Versatile full-vectorial finite element model for crossed gratings,” Opt. Lett.34, 2216–2218 (2009). [CrossRef] [PubMed]
A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements (Springer, New York, 2004).
L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A14, 2758–2767, (1997). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper