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

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.


Introduction
Since their introduction by Bérenger in [1] 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,3]) allow to consider PMLs as a complex change of co-ordinates, which implies equivalent material properties ( [4,5]). 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,7]), 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,9,10]). 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.

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 0 = ω/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 : where the coefficients ε xx , ε aa ,...µ zz are (possibly) complex valued functions of x and y, and where ε a (resp. µ 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 θ 0 : k + = α + x + β + y = k + (sin θ 0 x − cosθ 0 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) : where A 0 e (resp. A 0 m ) 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 : such that the diffracted field E d := E − E 0 and H d := H − H 0 satisfies an Outgoing Wave Condition (OWC, [11]) 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 µ : the functions e and h are solution of similar differential equations : with for the TE case, 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 ε + µ + , where k 0 := ω/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 ε − µ − • 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.

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]). 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]) : 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 : 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 C 3 , which co-ordinate change leads to equivalent permittivity and permeability tensors. We detail here the different co-ordinates used in this paper.
• (x s , y s , z s ) are the complex stretched co-ordinates. A suitable subspace Γ ⊂ C 3 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 = Re(x s ), y = Re(y s ), z = Re(z s )).
• (x c , y c , z c ) are three real co-ordinates corresponding to a real valued parametrization of Γ ⊂ C 3 .
In this paper, we use rectangular PMLs ( [12]) 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 : The expression of the equivalent permittivity and permeability tensors are thus : 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 : (11) where ζ ± = ζ ′ ,± + iζ ′′ ,± are complex constants with ζ ′′ ,± > 0.
In that case, the complex valued function y(y c ) defined by Eq. (9) is explicitly given by : 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 : where Eventually, denoting β ± and accounting for an OWC, we obtain the expressions of the different diffraction orders : 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 : The non oscillating part of this function is given by : For this, two pitfalls must be avoided : 1. 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.
2. 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 : and we take the larger value among the l − n : We set the height of the bottom PML regionĥ − = 10l − .

The FEM formulation
Our FEM formulation is that described in [13,14]. 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 1 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]). In the numerical examples in the sequel, the maximum element size is set to λ 0 /(N m Re(ε 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. 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 0 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 h − , 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 λ 0 = 10µm and of angle of incidence θ 0 = 10°in the air (ε + = 1). All materials are non magnetic (µ r = 1) and the periodicity of the grating is d = 4µm. We setĥ − = 10l − 0 and ζ ′ ,− = 1.

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 → β − 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 : so that the application β − defined by β − (y d ) 2 = k 2 0 ε − − α(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 n d such that around this depth the apparent wavelength corresponding to the order in play is contracted to a value close to λ 0 . 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 : The contraction y(y d ) is thus given by : The function y(y d ) has two poles, denoted y ⋆ d,± = d(± √ ε − −sin θ 0 ). When y ⋆ d,± = ±nλ 0 with n ∈ N ⋆ , β − (y ⋆ d,± ) = β − (±nλ 0 ) = β − ± = 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,± associated with a Wood's anomaly, denoting : In a second step, we seek the point y 0 d = n ⋆ λ 0 such that : To avoid the singular behaviour at y d = y ⋆ d,± , we continue the graph of the function y d (y) by a straight line tangent at is the so-called stretching coefficient. The final change of co-ordinate is then given by : (21) Figure 4 shows an example of this co-ordinate mapping. Eventually, the complex stretch s y used in Eq. (10) is given by : 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.

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.999y ⋆ d,+ ).Moreover, we set the length of the PMLĥ − = 1.1y ⋆ 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. 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). 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]). 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 0 + T −1 + T 0 + T +1 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. 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.

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.